A Physics-Informed Assembly of Feed-Forward Neural Network Engines to Predict Inelasticity in Cross-Linked Polymers

In solid mechanics, data-driven approaches are widely considered as the new paradigm that can overcome the classic problems of constitutive models such as limiting hypothesis, complexity, and accuracy. However, the implementation of machine-learned approaches in material modeling has been modest due to the high-dimensionality of the data space, the significant size of missing data, and limited convergence. This work proposes a framework to hire concepts from polymer science, statistical physics, and continuum mechanics to provide super-constrained machine-learning techniques of reduced-order to partly overcome the existing difficulties. Using a sequential order-reduction, we have simplified the 3D stress–strain tensor mapping problem into a limited number of super-constrained 1D mapping problems. Next, we introduce an assembly of multiple replicated neural network learning agents (L-agents) to systematically classify those mapping problems into a few categories, each of which were described by a distinct agent type. By capturing all loading modes through a simplified set of dispersed experimental data, the proposed hybrid assembly of L-agents provides a new generation of machine-learned approaches that simply outperform most constitutive laws in training speed, and accuracy even in complicated loading scenarios. Interestingly, the physics-based nature of the proposed model avoids the low interpretability of conventional machine-learned models.


Introduction
Exponential growth of the computational power over the last decade has enabled the first-generation of Machinelearned (ML) models to be used in computational mechanics [1] and polymer physics [2].Current ML models were often developed based on "black box" approaches, which beside low interpretability, require large volume of training data to prescribe certain behaviour.In solid mechanics, the stress-strain tensors are only partially observable in lowerdimensions, and thus providing data to feed black-box ML model is extremely challenging.Data-driven approaches in computational mechanics can be categorized into three groups (see review [3]) Model-free Distance-minimization Approaches .are developed to circumvent the need for analytical models by directly finding stress-strain points with least distance to the experimental data, which also satisfy the constraints set by kinematics and essential conservation laws.The approach was initially set for linear elastic materials [4], and later expanded to include hyperelastic materials [5] and different states of deformation.This approach suffers from a few issues; excessively high computational cost, strong sensitivity to data scattering, and higher-dimensionality problem induced by lack of data [6].
Non-linear Dimensionality Reduction Approaches .seek to build a constitutive manifold from experimental data to describe an accurate approximation of the strain energy in different states of deformation.The approach focuses on describing the constitutive behaviour through a set of predefined shape functions, such as β − spline [7], with 1 Corresponding author.constants derived by linearization strategy [8,9], or ML approach [10].Being mainly derived from WYPiWYG model [11], it focuses on solving the PDE equations obtained by shape functions, rather than fitting them.In elasticity, Manifold Learning is shown to be more efficient and more accurate than black-box ML models, and it has already been generalized to cover damage [12].The main problem with this approach is the large number of tests needed for validation, and its dependency on the assumption of constitutive manifolds with a particular functional structure [13].
Autonomous Approaches .incorporate ML models as surrogate functions to capture the high-dimensional and nonsmooth micro-scale behaviour of the material constituents, which has been shown to be a successful approach in Multiscale analysis [14].Several multi-scale analysis are proposed based on implementation of micro-scale ML model into the reduced-order FE simulations of macro-scale [15], a coupling which allows scalable utilization of ML surrogate models.However, the validity range of current ML models is extremely limited due to the unconstrained search space of the optimization variables, neglecting the problem's physics, difficulties in deriving parameter feasibility range, and lack of transition models to reduce the order of the problem.Implementing the reinforcement learning concept in meta-modeling of materials, a new class of ML models has been successfully developed based on cooperative/noncooperative game, where a pair of agents are trained to emulate certain performance through turn-based trial and errors [16].
Here, a cooperative multi-agent system A i j , i ∈ {1, n} , j ∈ {1, m} is proposed to describe different features in the material behaviour with n × m different learning agents.To reduce problem dimensionality, the 3D matrix is represented by m 1D directions, which allows us to replicate each agent m times to represent the 1D behaviour of the material.Each agent is then trained to emulate a certain behaviour of the material with the objective function being the error between the overall prediction of the system and the experimental data.Model fusion is used to integrate all agents back into a centralized system.

Physics-based Reduction
Helmholtz free energy (Ψ) was introduced as a summation of internal energy and entropy, which is functions of both deformation and temperature.Most constitutive models are developed to predict mechanical behavior of material under isothermal conditions, which are only function of deformation.Materials for which elastic deformations are reversible (e.g.entropy generation and internal dissipation is zero), and whose stress state can be derived from a potential, are called hyperelastic.Also, for homogeneous materials with uniform distribution of their internal constituents on the continuum scale and identical material properties in every material point, the position vector can be dropped as an argument.So, constitutive equations, for hyper-elastic materials, are derived directly from Clausius-Planck form of second law of thermodynamics which can be expressed using different work conjugate pairs of stress and deformation tensors (e.g.two-point tensors (F/deformation gradient:P/first-order Piola stress), material tensors (E/Lagrange strain:S/second-order Piola stress), spatial tensors (L/Hencky strain:τ/Kirchhoff stress)).In case of high polymer density, entropic energy is dominant part of free energy; however, accuracy fails when density of material reduces in porous materials or low-density gels as the contribution of energetic part increases.Strain energy function must accompany some conditions such as normalization, growth conditions, isotropy, objectivity, and polyconvexity, which guarantees the existence and uniqueness of the solution (Supplemental Material).Considering the challenging task of experimentally measuring the stress, strain fields within the matrix in the course of deformation, most modeling approaches are mainly based on the limited information available from the collective sample behaviour or Digital Image Correlation (DIC) reconstruction of 2-D strain fields.To address the challenge of significant missing data in high-dimensional data-driven approaches, we proposed a physics-driven order reduction approach by coupling the concept of Micro-sphere, Network Decomposition, Continuum Mechanics, and Polymer Physics to develop a sufficiently constraint machine-learned model that can predict the material behavior mainly based on macro-scale collective behaviour of sample.Fig. 1 demonstrates schematic of proposed model simplification idea.

Kinetics
The strain-displacement relations are obtained from a variation in the positions of a body, as defined by some coordinate system, both before and after deformation.The usual development of the strain-displacement relations considers a line element in a field.Consequently, the system's internal state is subject to the compatibility and equilibrium constraints which can be written as where f is applied force, B is a strain operator, u is displacement tensor and J = detF.

Continuum Mechanics
The response of a material subjected to thermal and mechanical loads is specified based on the relation of the stress (σ), internal energy (Ψ m ), heat generation (q), and entropy (η).Thus, the constitutive equations for a thermo-elastic material can be derived based on thermodynamics laws (see Truesdell et al. [17]), where (P) can be calculated by the partial derivative of the Helmholtz free energy with respect to (F) It should be noted that different conjugate pairs of stress-strain tensors ((i) two-point tensors, (ii) material tensors, and (iii) spatial tensors) can be used in Eq. 2. Hereafter, by enforcing polymer thermodynamics, the learning agents are first constrained to describe Ψ m (F) only.This condition automatically ensures the material objectivity, thermodynamic consistency, and polyconvexity requirement of the proposed constitutive model, as presented in the Supplemental Material.

Micro-Sphere
Polymer chains, in amorphous polymers, are arranged uniformly distributed within their polymer matrix.Due to non-specific spatial arrangement of polymer chains, researchers utilized micro-sphere concept for scale transition to model amorphous polymer network.In these models, unit-sphere is a well-known concept to reduce the order in material models by assuming a 3D matrix to be a homogeneous assembly of 1D elements distributed in all spatial directions (see Fig. 1).This concept assumes a similar distribution of elements in all directions in the virgin state, a condition which evolves based on the applied loading in different directions.Meanwhile, if coupled with strain amplification approach, it addresses the problem of inhomogeneity in micro-stretch distribution because stretch of the polymer chains between aggregates (micro-stretch) generally exceeds the stretch applied to the polymer matrix.To address this challenge, a concept of hydrodynamic reinforcement can be used, which, can be integrated into NN.This approach can transfer information from micro-structure behavior to the macroscopic behavior via homogenization over the unit-sphere.In micro-sphere approach, the average response of the material over the sphere can be numerically calculated by n integration directions [d i ] i=1...n that are weighted by factors [w i ] i=1...n [18].Hence, strain energy function over the sphere can be written as Step 1 Step 2 Step 3 Step 4 Here, all directions have the same strain energy function, which make this model proper for isotropic materials which their invariants remain constant during rotation.
The micro-sphere concept addresses the challenges of struggling with 3D tensors by changing the problem to 1D sub-networks.Also, because of the directional dependency of the model, it is more applicable in predicting the onset of fracture, deterioration, and propagation of cascading failure, which is extremely sensitive in history-dependent materials.Besides, we can apply complex loading by employing this model.

Network Decomposition
In the previous sections, the prediction of 3D energy function of material is simplified to a 1D problem.For further simplification, the 1D networks is decomposed to several sub-networks, where each of them captures a nonlinear feature in the response of material as a simple NN.Thus, strain energy in each direction is the result of summation of sub-networks' energy in given direction (Fig. 1).These sub-networks act in parallel ( . So, the total strain energy of micro-sphere can be written as where ψ j d i indicates strain energy of sub-network j in direction d i , and N s is the number of sub-networks. Energy of each sub-network is the output of NN related to that sub-network, ψ d i j = NN j (W j , S j ), which W j denotes weight matrix, and S j is vector of inputs designed for j th sub-network based on the problem.The input vector can be different features of the material including stretch, time, temperature, and etc. as input.Consequently, based on Eqs. 2 and 4, the first Piola-Kirchhoff stress tensor P can be written as where p denotes multiplier Lagrange to guarantee incompressibility.To train the model, the cost function should be minimized subjected to set of physical constraints described in section above.The cost function can be written as which P n shows the experimental data of stress with index n in data set.
Training Agents: To train 1D behavior of sub-networks, NN-based learning agents which demonstrate energy of sub-networks in an arbitrary direction are designed based on several inputs and internal non-kinematic parameters.These inputs are chosen based on behavior and complexity (e.g.full memory/recent memory) of materials A i j := ψ i j (E, M) where E and M represent the set of input and internal non-kinematic parameters.Also, for more complex mapping, we need to chose a greater number of layers and neurons; however, even one layer is sufficient in most of the cases.Note that training agents should satisfy normalization, growth conditions, isotropy, objectivity, as well as polyconvexity.
Material With Full/Recent Memory: For history-dependent materials, internal parameters of agents are designed based on type of materials' memory.For materials with recent memory like visco-elastic materials, internal memory of each sequence affects next sequence.However, in material with full memory such as rubbery materials with damage, only the maximum status of history affects next sequence.Training agents with full history is usually much more expensive since the response highly depends on a set of events.

Implementation to Rubber Inelasticity
To show the performance of devised model, inelastic and nonlinear behavior of rubber materials are investigated in the following.We need to define unknowns to solve the problem.The number of directions in a micro-sphere model shows the accuracy and efficiency of integral estimation over micro-sphere, which was investigated in this research [18].The greater number of directions is, the more accurate integral estimation will be.To show accuracy of the derived model, we choose the least amount of directions (21 directions).To design the number of sub-networks, we consider that sub-networks should be representative of first and second invariants to predict different states of deformation simultaneously [19].Besides, the inputs and internal parameters of learning agents are designed to capture the deformation of materials with full memory for rubbery materials.So, the inputs and internal parameters of full memory which we show through • max parameters will be where C = F T F is right Cauchy-Green deformation tensor.Final step is design of NN in learning agents.To have a simple agent trained fast and create complex mapping, we consider a multilayered NN.By trying different structures, we found that a NN with one hidden layer consists of four neurons with non-linear activation functions can capture mechanical behavior of this material well.Here, the cooperative multi-agent system will be A i j , i ∈ {1, 21} , j ∈ {1, 2}.Thus, loss function after the fusion of agents can be written as subjected to weights related to λ 1max and λ 2max 0 ; and weights related to λ 1 and λ 2 0 to satisfy thermodynamic consistency and polyconvexity respectively.Fig. 1 shows the schematic concept of the derived model.
Data-set Minimization.The selection of stress-stretch data set plays a vital role in the success of the model.This data set for training should cover the range for which this model will be used.Due to the contribution of compressed directions in uni-axial tension, we can say that the second NN in a sub-network assembly is negligible for uni-axial tension training.The physical explanation for this phenomenon is that with the same maximum stretch for uni-axial and bi-axial tension, the value of sub-networks' compression for different directions in bi-axial tension is much larger than uni-axial.Therefore, with the training of uni-axial data, we cannot capture bi-axial behavior correctly because NN cannot learn the difference of compression between these modes.To train the model correctly, we can consider three options.First option is training with bi-axial data.Fig. 2.a and Fig. 2.b shows the difference in performance between training with bi-axial and uni-axial data.The second option is that we can train with uni-axial data that its maximum stretch is larger than maximum stretch in bi-axial tension because, with this, we can compensate for the difference of sub-network compression between uni-axial and bi-axial tension.Fig. 2.c shows the result of this phenomenon.The last option is, training the model with uni-axial compression and extension to predict the behavior of shear and bi-axial tension Fig. 2.d.Table 1 shows the range of prediction related to data set training type.[23] To show the model can predict non-linear and inelastic behavior of rubber for multi states of deformation, several experimental data sets are employed from literature.The model is trained with bi-axial tension and predict uni-axial tension and shear (see Fig. 4 and Fig. 5).For observing performance of the model with more results please see Supplemental Material.

Conclusion
Here, we introduced a new physic-informed data-driven constitutive model for cross-linked polymers by embedding NN sub-networks into micro-mechanical platform.This is a less data-dependent, fast trained, low dimension, and interpretable model such that micro-mechanical behavior of sub-networks in each direction is obtained, including all interactions, directly from macroscopic experimental data set.We employed this approach to indicate performance and excellent accuracy of the model by using several experimental data sets such as uni-axial extension and compression, pure shear, shear compression, and bi-axial extension for rubbers.The model demonstrates significant performance in prediction of both deformation state dependency and inelasticity, such as Mullins effect and permanent set.To the best of our knowledge, this is the first time such an achievement has been met.Also, physical constraints, such as polyconvexity, frame independency, and thermodynamic consistency, have been satisfied.

Figure 1 :
Figure 1: Schematic of proposed model for a micro-sphere with two sub-networks

Figure 3 :
Figure3: Model training and prediction of cyclic uni-axial tension with step-wise increasing of amplitude[23]

Table 1 :
Prediction domain for train till stretch χ