Generalizing the Gurson model using symbolic regression and transfer learning to relax inherent assumptions

To generate material models with fewer limiting assumptions while maintaining closed-form, interpretable solutions, we propose using genetic programming based symbolic regression (GPSR), a machine learning (ML) approach that describes data using free-form symbolic expressions. To maximize interpretability, we start from an analytical, derived material model, the Gurson model for porous ductile metals, and systematically relax inherent assumptions made in its derivation to understand each assumption’s contribution to the GPSR model forms. We incorporate transfer learning methods into the GPSR training process to increase GPSR efficiency and generate models that abide by known mechanics of the system. The results show that regularizing the GPSR fitness function is critical for generating physically valid models and illustrate how GPSR allows a high level of interpretability compared with other ML approaches. The method of systematic assumption relaxation allows the generation of models that address limiting assumptions found in the Gurson model, and the symbolic forms allow conjecture of decreased material strength due to void interaction and non-symmetric void shapes.


Introduction
Constitutive models are often developed using a combination of experimental observation, mechanics, and physics-based principles.In these derivations, cost and fidelity balances can be achieved by idealizing otherwise complex microstructures.These idealizations can limit model accuracy and generality to smaller subsets of materials or loading conditions.A prevalent example of this outcome is the Gurson model: a micromechanics-derived formulation for ductile rupture in which various assumptions are introduced with idealized microstructures to enable analytical derivation [1,2].Extensions to reduce the limitations imposed by such idealizations can be made by adding adjustable parameters that are optimized to fit data from either experiment or numerical simulation [3].These models gain accuracy with respect to real-world materials but add parameters that must be optimized for each new material to which the model is applied.In addition, these added parameters provide little physical insight into the mechanics of the material, limiting model interpretability.
An alternative to relying on idealizations that permit closed-form derivation is to use machine learning (ML) to train models using experimental or simulation data.There are many examples of studies using artificial neural networks (ANNs) for constitutive modeling [4][5][6].However, these ML approaches lack inherent interpretability, which has been a primary benefit of the more traditional material model development methods.Physics-informed ML (PIML) is often used to generate more explainable and physically valid material models.PIML relies on incorporating known partial differential equations as well as numerical constraints and limits into the ANN loss function [7].Specifically, in the case of constitutive modeling, incorporating thermodynamic principles into an ANN can lead to more efficient training and thermodynamically consistent model forms [8].While these approaches allow the generation of models that better abide by the known physics of the system, they are still limited in their interpretability due to the black-box nature of ANNs.
Interpretability is important in material modeling for improving understanding and trust in the constitutive models these algorithms generate.By generating inherently interpretable models, i.e. in the form of mathematical expressions, engineers can directly gain an understanding of the material by seeing how features are incorporated into the model form.This direct understanding also helps improve trust in these models, as they can be related to known physical principles.An additional advantage is that models in these forms can be readily integrated into existing engineering workflows such as finite element (FE) analysis.Finally, increased model interpretability can help extend existing knowledge and allows the incorporation of human intervention into the model training process.
An inherently interpretable ML approach, genetic programming based symbolic regression (GPSR), generates free-form symbolic expressions from training data [9].Previous work utilizing GPSR for constitutive modeling has shown promise in generating accurate constitutive models; however, interpretation of these models for complex microstructures has been limited.Multiple studies have utilized symbolic regression (SR) to replace tunable parameters within constitutive models with microstructural feature-dependent equations [10][11][12].Versino et al [13] used SR to generate a flow stress model within the von Mises plasticity framework based on experimental stress-strain data for copper.In Bomarito et al [14], GPSR was trained on FE simulation data of a porous ductile metal to generate the equation for the yield surface of the material.The resulting model was compared to the Gurson model for ductile porous metals [1,2].Interpretability of the final output model was limited since multiple idealizations of the microstructure had been made.
Starting from the Gurson model, we seek to maintain the interpretability of GPSR material models upon consideration of increasingly complex deformation mechanisms by systematically relaxing Gurson's assumptions.By relaxing assumptions in this way, we can directly relate additional terms in the output GPSR models to the assumptions and microstructure idealizations.Several transfer learning (TL) methods are used to complete the assumption relaxation steps and their results are assessed.Any of the presented TL methods can be included within an automated process, but some methods demonstrate advantages which are also discussed.In section 1.1, we provide a background on ductile failure modeling of porous materials, which will be the focus of this work, followed by the assumptions present in the Gurson model.In section 1.2, we will discuss the many approaches taken to improve the model regarding these limiting assumptions.

Modeling of porous ductile materials and the Gurson model
Ductile fracture in metals typically occurs by nucleation, growth, and coalescence of microscopic voids within the material [15].Modeling the growth of these voids and how they affect material performance is at the center of ductile failure modeling.The first micro-mechanical studies on the phenomenon were completed in 1968 when McClintock studied the evolution of cylindrical holes in a plastically deformed material [16].Shortly after, a similar study was conducted with spherical holes [17].McClintock [16] and Rice and Tracey [17] noted the dependence of void growth on the hydrostatic component of stress, a critical mechanical characteristic of porous metals.Soon after this early work, Gurson completed micro-mechanical studies and averaging techniques from Bishop [18] to develop an upper bound yield criterion for porous ductile materials [1,2].This work proved to be foundational on the subject of ductile failure and has since been used extensively in various applications [15].
The Gurson model for porous ductile materials was formulated using a hollow sphere and a hollow cylinder subjected to a homogeneous boundary strain rate.The materials were idealized as rigid-perfectly plastic von Mises materials.Gurson used an upper bound approach to derive the formula for the yield surface, prescribing a macroscopic flow field and solving a dissipation integral to find the minimum dissipation for the model [1,2].This upper bound ensures that the predicted yield strength is always greater than what occurs in the material.The void volume fraction (VVF) was used to describe the level of porosity, defined by the volume of the voids normalized by the corresponding volume of the material.Equation ( 1) is the resulting yield surface equation: where σ vm is the von Mises stress, σ y is the matrix material yield strength, σ h is the hydrostatic stress, and f is the VVF of the material [1,2].A characteristic of this yield surface is that the equation simplifies to that of the von Mises yield surface when VVF → 0. This behavior is expected because the derivation assumes a rigid-perfectly plastic von Mises matrix material.
Along with the yield surface, the evolution of the VVF is found through the conservation of mass, resulting in the void evolution equation: where εp is the macroscopic plastic strain rate and bold notation is used to reflect variables that are tensors.Equations ( 1) and ( 2) define the Gurson model of yield and void growth in ductile metals, respectively.The Gurson model has proven to be accurate and tunable for modeling porous metals in a wide variety of problems [15].However, microstructural idealizations and assumptions made in the Gurson model formulation limit accuracy and generalizability to diverse stress states such as those with low triaxiality.These limiting assumptions are: 1.The theoretical model is based on a single void, which does not account for interactions that can occur between neighboring voids in a material.2. Voids within the material are idealized as spherical and must remain spherical throughout deformation, which is especially unrealistic under deviatoric deformation gradients.3. The model does not consider randomly placed voids within the material, precluding the effect of void clusters and their influence on behavior variability.4. The material is idealized as perfectly plastic, meaning no hardening is present, which is not representative of real-world material behavior.5.The model does not account for void coalescence within the material, which limits its ability to predict material behavior under large deformations.6.The model does not account for void nucleation, which limits its predictions for void growth in materials with small initial VVFs.
In order to generate improved, interpretable models that eliminate these microstructural idealizations, we will utilize GPSR trained on FE simulations where these limiting assumptions are relaxed or removed one at a time.Thus, newly generated models will be directly correlated to the assumptions present in the material, improving model explainability.In this work, we will focus on elimination of the first two assumptions.

Review of Gurson model improvements
An important series of improvements to the Gurson model led to what is known as the Gurson-Tvergaard-Needleman (GTN) model [3,[19][20][21].This modified Gurson model added tunable parameters within the yield surface function as shown in equation (3), extending the Gurson model to a wider variety of materials and loading cases.The GTN model also considered void growth due to nucleation.Despite increasing model generality, the additional parameters did not help explicitly address the assumptions listed in section 1.1 To understand the limiting effect of assumptions 1 and 3, Fritzen et al [22] developed an FE model of a representative volume element (RVE) with randomly distributed voids and compared simulated results to that of the Gurson model.They concluded that the Gurson model over-predicted yield compared to the FE model for larger values of VVF (values over 1%), providing direct evidence that the Gurson model is an upper bound and demonstrating the effect of assumptions 1 and 3. Bilger et al [23] conducted a similar study on microstructural models containing random void locations utilizing fast Fourier transform (FFT) numerical simulations.They were able to show the behavior of shear bands for microstructures with and without void clusters and compare the resulting yield surfaces to the Gurson model.Similar to Fritzen et al [22], the comparison showed that the FFT models predicted yield at lower stress values than Gurson [23].
The limitation of assumption 2 on Gurson model predictions is most prominent for low triaxiality load conditions (high von Mises stress and low hydrostatic stress, occurring near pure shear) [24].These load cases exhibit void shape change through deformation, not accounted for by the Gurson model due to assumption 2. Nahshon and Hutchinson [25] attempted to resolve this shortcoming by altering equation (2) to account for the deviatoric strain.Inclusion of the Lode parameter, a function of the third stress invariant, helps indicate the character of the stress state and leads to better void growth predictions in regions of low stress triaxiality [25].Danas and Castañeda [26] and Daehli et al [27] developed similar methods to incorporate the Lode parameter into the Gurson model, fixing predictions made by Nahshon and Hutchinson [25] for material ductility given both high and low stress triaxialities.
Additional research efforts have been performed to address the limitation of assumption 4 by attempting to modify the Gurson model to account for hardening in the matrix material.Leblond et al [28] analyzed both the isotropic and kinematic hardening case, performing an analysis similar to Gurson's with a hollow sphere of some hardenable material.As a result, the σ y term inside the square and cosh terms of equation ( 1) was replaced by two new variables: σ 1 and σ 2 .These new variables are the effective yield stress under deviatoric and hydrostatic loading conditions, respectively.By modifying the Gurson model, predictions of yielding and void growth for purely deviatoric and hydrostatic loadings were improved when compared to FE simulations [28].Another implementation of kinematic hardening was introduced by Mear and Hutchinson [29], who implemented a macroscopic back-stress into Gurson's formulation.The addition of a back-stress resulted in a significantly more complex evolution equation than Leblond [28].More recently, work was conducted to create a Gurson-type model for both isotropic and kinematic hardening and was created by analyzing a hollow sphere in discretized layers to help account for the heterogeneity of hardening from the center to the outside of the sphere [30].
Work on incorporating void coalescence into the Gurson model to account for violations of assumption 5 started with the GTN model, where void coalescence was included as an acceleration in void growth past a particular critical value of VVF, f c as shown in equation ( 4) [3,[19][20][21]] The critical value of VVF could then be tuned to match experimental data for any porous material.Additional work has been done to estimate this critical VVF value based on analyses of shear localization in Gurson-like materials [31].Finally, FE model calculations of RVEs have been used to infer these parameters to some promising levels of success [32].These improvements highlight how each assumption made in the Gurson model may lead to deficiencies in its predictive capabilities.Many studies also present ways of mitigating these issues.Using an interpretable ML approach, we will present a new method to generate models with fewer assumptions.These previous studies will guide the GPSR training process and lend insight into the mechanics of produced model forms.The gradual relaxation of assumptions will allow a direct connection from the model form to the relaxed assumption to increase physical insight and interpretability.

Methods
To improve the interpretability and explainability of material models using GPSR, we relax the assumptions in the Gurson model, one at a time.This systematic relaxation requires controlling the assumptions in the training data, so FE simulations were used instead of experimental data where controlling these assumptions would be impossible.We then train GPSR models on each dataset, each sequentially relaxing an assumption.We start with dataset 1, which relaxes assumption 1, then dataset 2 relaxes assumptions 1 and 2, and continue until all assumptions are relaxed.Here, we will present results for the first two assumption relaxations listed in section 1.1.This assumption relaxation process also lends itself to incorporation of prior knowledge through TL to maximize training efficiency and model interpretability.This process of relaxing assumptions and producing GPSR models is illustrated in figure 1.

Training data generation
An FE model is generated for each assumption relaxation step to produce corresponding training data, table 1, using Abaqus/Standard [33]. 4 A unit cube of material with a single void embedded in the matrix is used for the first and second assumption relaxation.To enforce assumptions 2-6, the nodes surrounding the void are constrained so that the void remains spherical throughout deformation.These constraints on the shape of the void are referred to as the self-similar constraint in table 1.To relax assumption 2, the self-similar constraint on the void is removed.The material of the RVE is consistent with the Gurson model formulation, with a matrix of perfectly plastic von Mises material.The material properties required for this model are a Young's modulus, E = 200 GPa, a Poisson's ratio, ν = 0.3, and a yield stress, σ y = 100 MPa.These material properties are chosen to resemble a material that fails via mechanisms consistent with those of ductile fracture.The two FE models used to generate the datasets for the assumption 1 and 2 relaxations, respectively, are described in table 1.The FE models for the assumption 1 and 2 relaxations are meshed using quadratic tetrahedral elements, and convergence studies were conducted to ensure the mesh was sufficiently converged.
Meshes with a relative error less than 0.05% for VVF and 0.5% for homogenized von Mises stress are used for training simulations.An example of the RVE modeled in Abaqus is shown in figure 2(c).Contact between surfaces of the void is not considered in the FE model since the applied deformations for assumptions 1 and 2 are small enough that no contact occurs; however, addition of contact will be required once further relaxations, like void coalescence, is modeled.Load cases are generated using a space-filling technique outlined by Zamora [34], designed to distribute datapoints across the chosen stress space equally.This technique results in a variety of applied boundary conditions (displacements in the x, y, and z directions) to produce deformations spanning shear-dominated, uniaxial, and near-hydrostatic loading.Here, we use 200 triaxial load cases generated using Zamora's [34] technique that are applied via periodic boundary conditions (PBCs) to all six sides of the RVE.The PBCs work by applying displacements at four control nodes, figure 2(c).Constraint equations between the control nodes on opposing faces and the rest of the corresponding face nodes ensure opposing faces of the RVE are under the same imposed strain [35].The void size is changed to evenly distribute the datapoints between an initial VVF (f in figure 2(a)) of range 0%-10%, and the VVF evolves according to the applied deformations.The amount of influence that the assumptions 1 and 2 relaxations have on the data depends on the size of the void and the deformation gradient being applied.For example, as opposed to an assumed infinite medium surrounding a void, PBCs permit interactions among voids, the effect of which increases with VVF.The deformation gradients of each load case are discretized evenly over 145 load steps to generate training data dependent on plastic strain and state variable evolution.The load step discretization was selected based on a convergence study similar to the one used to determine mesh size.Variables of interest are homogenized from the RVE using a volume average of values taken from each element and stress-like variables are normalized by the matrix yield strength, σ y .The variables of interest are those found in the Gurson model shown in equation ( 1) and the equivalent plastic strain to track when the model has entered the plastic regime.The resulting dataset is where the overbar on ϵ p is used to signify that it is the equivalent plastic strain.Each variable in equation ( 5) is treated as a vector of values across the 145 load steps.Since the desired GPSR Points are colored according to their equivalent plastic strain to indicate points at the onset of plasticity (low equivalent plastic strain) and end of the load case (high equivalent plastic strain).Finally, it is observed how the yield surface constricts at larger values of f, indicating the expected decrease in overall strength as more and more material is removed.

GPSR algorithm
Once training data is generated using the FE model, we use GPSR to find an equation for the yield surface model.Potential models generated by GPSR are written as symbolic equations are represented as acyclic graphs and are implemented in the code as lists of operations, figure 3. The equation shown in figure 3 is only used to illustrate how equations are represented in GPSR and does not relate to the dataset described in the previous section.Once the models are defined, genetic programming (GP) evolves models that best fit the data.GP evolves new models using two methods: crossover and mutation.Crossover combines part of an operation list of a model with part of an operation list of another model.Mutation either (1) deletes an operation, (2) changes the operation and parameters, or (3) only changes the parameters of an operation in an operation list of a model.The probability that crossover or mutation will occur, p cross and p mut , respectively, at each evolutionary step is defined by a hyperparameter.Crossover and mutation evolve models in the population, and models with the best-defined  fitness are selected for subsequent generations.For a more extensive description of GPSR, see Schmidt and Lipson [9].Models evolve with undefined constants, such as c0 shown in the acyclic graph in figure 3.These coefficients are then calibrated to the training data using the Levenberg-Marquardt optimization algorithm [36].GPSR models are trained using the open-source NASA code, Bingo [37].Bingo uses the coevolution of fitness predictors and island population division to speed up GPSR training [37].Coevolution of fitness predictors evolves a representative subset of data for fitness evaluations to reduce the computation cost of fitness evaluations [38].Island population division separates the population of models across multiple computer processors, referred to as islands, allowing them to evolve in parallel.Individuals from separate islands are randomly migrated at a user-specified migration rate.The hyper-parameters used for this study were the same as those used by Bomarito et al [14] and are shown in table 2.

Fitness function.
GPSR selects equations based on minimization of their fitness with the data.Defining fitness poses a challenge when the evolved equations of interest are implicit, i.e. they do not have a dependent variable.Implicit equations are of the form F(X) = c, where c is a constant.Using a traditional explicit fitness metric, such as mean squared error, would lead to arbitrary solutions, e.g.x 0 − x 0 + c = c.To avoid arbitrary solutions, Bomarito et al [14] devised a fitness metric for implicit equations based on the work of Schmidt and Lipson [39] that compares the partial derivatives of a model with those of the data.The metric provided by Bomarito et al [14] was used to measure model fitness in this study: ∂F ∂xj where E is the model fitness, F is the model, x j is an individual variable of dataset X, N is the number of data points, P is the number of parameters, and t is the pseudo-time.The metric is based on the Prager consistency condition, which states that dF/dt = 0 for models where F(X) = c.The chain rule is used to calculate dF/dt by taking the partial derivative of the proposed model with respect to each parameter and multiplying that by the change in that parameter over time.Therefore, dF/dt should be minimized to zero to abide by the Prager consistency condition.Although the training data used will always abide by the Prager consistency condition, not all proposed GPSR models will abide.The absolute value in the denominator of equation ( 6) is used to avoid division by zero for proposed GPSR models that satisfy the Prager consistency condition, but at the same time avoid models that arbitrarily satisfy the condition, e.g.
If any values of not a number (NaN) are returned because of division by zero, the NaN values are removed from the fitness calculation.If more than 10% of the N calculations are NaNs, the fitness is returned as infinity [14].Using this definition of fitness avoids arbitrary solutions that would arise using a traditional explicit fitness metric, so this formulation is crucial for using GPSR with our training data.

Methods of prior knowledge transfer.
To improve model interpretability, we seek to implement TL to include known terms from existing models or knowledge in the GPSR training process.Zamora [34] established that the Gurson model can be found using GPSR, so we look to extend that work using TL to generate models with fewer assumptions than Gurson.To do this, two methods are utilized, referred to here as seeding and boosting.Seeding in GPSR originates from a study conducted by Schmidt and Lipson [40], wherein terms from a known model are embedded within the initial GPSR population of models.By contrast, in conventional GPSR, initial populations of models are randomly generated.Therefore, in addition to improved interpretability, seeding aims to accelerate evolution through a more informed initial population.Here, the Gurson model, equation (1), is used as the seed for the assumption 1 relaxation since the first TL step is building from that existing knowledge.The assumption 2 relaxation will then utilize seeds from the GPSR model generated using the assumption 1 relaxation dataset.Per suggestions from Schmidt and Lipson [40], 10% of the initial GPSR population is replaced with seeds.
Boosting is another method to instill prior information into the GPSR training process by weighting data points based on a previously generated model's fitness.The algorithm for this method of boosting is shown in algorithm 1, where F is the proposed GPSR model, A is the assumption being relaxed, W is the weights applied to the training data, E is the implicit fitness, equation (6), and [Z 0 , Z 1 ] is the range between which the weights, W, are normalized.The subscript b indicates the bth boosting step corresponding to the bth assumption relaxation, section 1.1.

Algorithm 1. Boosting with GPSR.
Given input data X and proposed model F b (X) for assumption relaxation A b : 1. Calculate weights using previous model .
First, we evaluate ∂F b−1 /∂t using part of equation ( 6) and assign it to the weight vector, W b,i .Second, the weights are normalized between Z 0 = 0.5 and Z 1 = 1.This normalization was conducted to ensure that data points were weighted evenly.If Z 0 is set to 0, points especially poorly represented by previous models would dominate the training and lead to solutions that did not represent the training data as a whole.Finally, the fitness for the proposed model is calculated as the weights multiplied by the implicit fitness metric, emphasizing the points that are fit worst by the previous models.Formulating the fitness this way ensures that newly formed models target these data points and compensate for previous model weaknesses.For the first two assumption relaxations presented in this study, the Gurson model was used to calculate the weights in step 1: where G(x) given in equation (1).

Fitness regularization.
After some preliminary testing, it was clear that many of the models produced by the seeding and boosting methods did not satisfy known constraints.For example, as noted in section 1.1 for the Gurson model, when VVF → 0, the model should reduce to the von Mises yield surface: σ 2 vm = 1, herein referred to as the von Mises constraint.The von Mises constraint should be satisfied for the first two assumption relaxations since the modeled material is still assumed to be a perfectly plastic von Mises material.To promote GPSR models that satisfy the von Mises constraint, a regularization constant, λ, is applied to equation (6) for models that match the von Mises surface at f = 0.The value of λ and its effect on the convergence and produced models is investigated.The complete algorithm is shown in algorithm 2.
1. Propose a model F(X) 2. Calculate model fitness E(F(X)) based on equation ( 6

Training comparisons
The three methods described (seeding, boosting, and fitness regularization) plus a conventional implementation of GPSR are all tested on the assumption 1 dataset, table 3. The best performing algorithms are then tested on the assumption 2 dataset.In addition, multiple penalty terms, λ, are tested for the fitness regularization algorithm, and the effect of different supplied operators is tested on the conventional GPSR algorithm.The four main training characteristics of interest are the model convergence rates, fitness, complexity, and interpretability.

Results and discussion
This work aims to develop a framework that automatically generates new plasticity models that improve interpretability for complex microstructures.By systematically relaxing assumptions associated with the Gurson model so that additional terms in the resulting model can be related to the corresponding assumption, new plasticity models can be generated.Additionally, prior knowledge is incorporated into the GPSR training process via TL to investigate how this affects model interpretability and rate of convergence of the GPSR algorithm.To assess these different methods, we look at the rate of convergence of the GPSR algorithms and the GPSR model forms compared to the Gurson model.All algorithms shown in table 3 were run ten different times until they reached 10 000 generations on the assumption 1 dataset.The best-performing algorithms were then used on the assumption 2 dataset.

Assumption 1 models
The first characteristic of interest of the proposed algorithms is their effect on the convergence rate of the GPSR algorithm.Each algorithm is run ten times, each time for 10 000 generations, and the best model fitness is taken every 100 generations and averaged across the ten runs.The resulting convergence plot is shown in figure 4. The initial fitness value of the seeding algorithm starts at the lowest or best fitness of all the algorithms.A low initial fitness illustrates the desired effect of using the Gurson model in the initial population, as it helps the GPSR algorithm achieve better model fitness early in model evolution and improves the overall rate of convergence.After 10 000 generations, however, the effects of seeding are minimal, as the boosting and conventional algorithms can achieve very similar values of fitness.Another interesting comparison is the conventional algorithm with and without hyperbolic cosine as an operator.When the hyperbolic cosine term is included in the operator set, the conventional GPSR algorithm can achieve better fitness since it enables the GPSR algorithm to more easily find equations that fit the training data.Without the hyperbolic cosine operator, the algorithm is forced to construct an equivalent hyperbolic cosine term using exponentials, cosh(x) = (e x + e −x )/2, which it is rarely able to do in less than 10 000 generations, leading to worse performance.The importance of setting the correct penalty parameter, λ, in the fitness regularization algorithm, algorithm 2, is illustrated in figure 4. With a penalty of 20%, the algorithm stagnates at a worse fitness than even the conventional GPSR algorithm without cosh.However, with a penalty of 5%, the algorithm can achieve similar levels of fitness to the conventional, seeding, and boosting algorithms.When the penalty is too large, the GPSR algorithm starts to overly focus on the von Mises constraint and fails to fit the dataset as a whole.Finally, the boosting algorithm appears to have minimal effect on the rate of convergence of the GPSR algorithm compared with the other algorithms tested, as shown in figure 4. The fitness difference between the best and worst models appears insignificant when visualized in stress space.This means that the algorithms produce models with nearly equivalent fitness values, but differences lie in their convergence behaviors and the forms of the models they generate, which is discussed below.
To understand the improvement in model accuracy at this assumption relaxation step, we compare a selection of models to the Gurson model.In a single GPSR run, we are left with up to 26 potential models representing the best-fitting model at each complexity level, where complexity is the number of operations in an operation list of a model, figure 3. The number of potential models depends on the number of equations at unique complexities that are present in the final equation population.In our trials, the maximum number of unique complexities was 26.These 26 models multiplied by the ten trials conducted for each of the six algorithms results in up to 1560 potential models across all algorithms and trials.One common approach to model selection is to choose the model that best fits the data.However, our approach seeks to choose models that fit the data, are interpretable, and lend new information about the material we are modeling.To select the best models following these criteria, we first identify models from each algorithm that abide by the von Mises constraint defined in section 2.2.3, ensuring that the models we choose abide by a known behavior of the material, as f → 0.
Enforcing this constraint left 99 models: 76 from the fitness regularization algorithm, 16 from the conventional algorithm, 5 from the boosting algorithm, and 2 from the seeding algorithm.Not surprisingly, the fitness regularization algorithm was the most crucial in generating physically valid models.Many of the 99 models have the same complexity, so to further reduce the number of models, only the models with the lowest fitness at each complexity were kept, resulting in 19 models whose Pareto front is illustrated in figure 5.In figure 5, the most complex models are not necessarily the best fit since this Pareto front of models was taken from different algorithms and GPSR trials.Out of these 19 down-selected models, there are a few different approaches to select a final model to analyze and compare with the Gurson model.First, we can select the model with the lowest possible fitness, which is shown in equation ( 7) and plotted in figure 5 with green triangles as the 'Best-Fit' model that came from the fitness regularization algorithm Although this equation satisfies the von Mises constraint, it is difficult to interpret as it has a cosh 2 term and hyperbolic cosine terms being multiplied.Equation ( 7) also lacks the f 2 term in the Gurson model, equation (1).Many of the best-fitting models have similar issues, containing complex terms that limit our ability to explain the effects on the yield surface.These terms are due to overfitting within the GPSR algorithm and demonstrate that picking the best-fitting model is not always the best in terms of complexity and interpretability.Another approach that generally improves interpretability and reduces overfitting is using Pareto front analysis to select a model that balances fitness and complexity, i.e. near a point where added complexity has insignificant reductions of fitness.This model is shown in equation (8) and shown in figure 5 with yellow squares as the 'Elbow' model that came from the fitness regularization algorithm This equation has a complexity of 11, is less complex than equation (1), and achieves a better fitness value.The only difference between this model and Gurson is removing the f 2 term and adjusting coefficients in the model form.This approach in fitting coefficients is consistent with other work, such as the GTN model [3], equation ( 3), but provides minimal insight into the effects of the specific assumption being relaxed.The final approach we investigated was to leverage GPSR as a feature engineering tool to generate an engineered model using the most common terms in the final down-selected Pareto front.Terms, as defined here, are portions of the GPSR model below an addition or subtraction operator.For example, using the acyclic graph shown in figure 3, the two terms would be x0 and c0 sin(x1).
We first record the most common terms present in the models of the down-selected Pareto front, figure 5.After analyzing all of the models, the most used terms are added together.The choice to add terms instead of using an alternate operation is chosen based on our definition of a term presented in the previous paragraph.Each term is then assigned an unknown coefficient that is later optimized.This algorithm is illustrated in algorithm 3. Algorithm 3. Generate engineered equation.
Initialize T[t j ] = 1 end if end for end for Initialize engineered equation E = 0. for term j, T j , in best n terms do E = E + C j T j end for Optimize E to solve for unknown constants, C, using implicit fitness.
Of the 19 models in the down-selected Pareto front, five terms were repeated in more than five models, while the rest of the recorded terms only appeared once.These repeated terms were (σ vm /σ y ) 2 , f cosh( σ h σy ), f 2 , f, and (σ vm /σ y )f.These terms were combined as described in algorithm 3, which resulted in equation ( 9).This model is identified as the 'Engineered' model in figure 5 This engineered model has a fitness of 0.071, better than 14 of the 19 down-selected models, and is straightforward to interpret due to its similarity to equation (1).This method of constructing an engineered model also ensures repeatability to other problems since, for example, equation (7) would not evolve with every GPSR run.However, the terms used to build the engineered model are more likely to appear with every GPSR run.The engineered model is plotted and compared to the Gurson model in figure 6.Each model is plotted in the von Mises stress and hydrostatic stress plane for multiple values of f.In figure 6, the engineered model is shown to be able to fit the assumption 1 data better than the Gurson model, as expected, since void interaction is relaxed in the training data generation.The improved fit is also reflected in the model fitnesses, which evaluated to 0.096 for the Gurson model and 0.071 for the engineered model.Comparing the two model forms, the engineered model contains two additional terms, the (σ vm /σ y )f and f term.Algebraically, this constricts the yield surface at high von Mises stress and high f.We can see this performance improvement by analyzing the difference in implicit fitness for the assumption 1 dataset, figure 7. The largest improvements in implicit fitness occur in regions of high von Mises stress due to the addition of the two terms.Relating to the physics present in the material, we can see that allowing the voids to interact causes a reduction in strength in areas of high shear stresses and high f.The latter observation makes sense intuitively since larger f values would   be expected to increase void interaction and decrease yield strength.The relation to the von Mises stress is less obvious, but many works have noted the deficiencies of the Gurson model in areas of low stress triaxiality (high von Mises stress and low hydrostatic stress) [28].Many of these works cite this deficiency as a change in void shape in loading cases with high shear stresses.However, here we can see that this deficiency may also be due to void interaction since that was the only relaxed assumption in this dataset.

Assumption 2 models
Due to the success of the seeding algorithm on improving the initial GPSR population fitness in assumption 1, figure 4, we chose to investigate two different seeding approaches for the assumption 2 dataset.Here, the initial GPSR population is seeded with all of the models from the assumption 1 Pareto fronts that abide by the von Mises constraint, including the three equations introduced in section 3.1, herein referred to as the regularized seeds.The regularized seeds are compared to the best-fit seeds, consisting of the best-fitting models from the assumption 1 Pareto fronts, where the models do not have to abide by the von Mises constraint.The goal of testing these two methods is to assess how seeding different equations into the initial GPSR population affects the convergence of the seeding algorithm.The two seeding approaches are tested ten times, and the best fitness out of the population is recorded every 100 generations and averaged across the ten runs, figure 8.
From figure 8, initially, the best-fit seeds have a better fitness value than the regularized seeds, which is to be expected since many of the models were overfitting for the previous assumption 1 dataset.However, as evolution occurs, the regularized seed models can achieve a lower average fitness across the ten runs.There are two potential reasons for this result.One, the models seeded using the regularized seed approach were, on average, less complex than those seeded using the best-fit approach, which would allow further evolution to occur on those models without reaching the maximum equation complexity set at 30.Second, the models in the regularized seed approach all abided by the von Mises constraint, which may have helped overall model fitness as evolution occurred.The second reason is less likely since, in both seeding cases, the von Mises term was evolved out of the final model forms, leading to none of the models produced by either method abiding by the von Mises constraint.The omission of the von Mises term is caused by the implicit fitness function, equation (6), which only relies on partial derivatives and not explicit checks of model limits to evaluate fitness.Because of this, the fitness regularization algorithm was still used to generate the final models for analysis since that algorithm was the only one able to produce physically reasonable models that abide by the von Mises constraint.
Using the engineered model approach shown in algorithm 3, the final Pareto front of all algorithms was reduced, and the most prevalent terms were found.Combining these terms resulted in the model shown in equation ( 10) Equation ( 10) contains a few additional terms not present in equation ( 9), including an additional multiplier in front of the f cosh(σ h /σ y ) term of σ vm /σ y as well as an additional σ vm /σ y f term, but removes the (σ vm /σ y )f term.Compared to the Gurson model, these terms again aid in constricting the yield surface at higher values of von Mises stress and VVF, and the multiplier in front of the f cosh(σ h /σ y ) term has a similar effect, but only at a larger value of hydrostatic stress.Visualization of this yield surface compared to the Gurson model is shown in figure 9.As with assumption 1, the GPSR model for assumption 2 fits the training data better, especially in areas of low stress triaxiality (high von Mises stress and low hydrostatic stress) due to adding the terms shown in equation (10).When the assumption 1 model shown in equation ( 10) is applied to the assumption 2 dataset, the difference is difficult to visualize in stress space.However, if we compare the calculated implicit fitness on all of the datapoints using the two engineered equations for assumptions 1 and 2, equations ( 9) and (10), we can see how this model improves with respect to the additional relaxed assumption, figure 10.
The model in equation ( 10) achieves a fitness of 0.069, compared to a fitness of 0.088 from the assumption 1 model, showing a notable improvement with the addition of these two terms.The assumption 2 engineered model improves most in areas of higher hydrostatic loading (greater than around 1.8), as is indicated by the higher difference in fitness in those regions, figure 10.We would expect the largest differences between this assumption 2 dataset to the previous assumption 1 dataset in regions of low hydrostatic stress and high von Mises stress, where the most void shape change is expected to occur.Instead, we show that the difference in those regions is fairly low, with the largest differences occurring in the higher hydrostatic loading cases due to the addition of the σ vm /σ y f multiplier in front of the hyperbolic cosine term.The location of these differences is likely due to the range of VVF used in this study, with a maximum of 10%.These VVFs ultimately resulted in a minimal difference between the self-similar void dataset and the previous assumption 1 dataset, whereas an even larger simulated VVF, though unrealistic for most metals, would potentially have lead to a larger difference in these assumptions.Thus, little information was gleaned from the assumption 2 model; however, both models showed notable improvement over the Gurson  model, particularly in areas of low stress triaxiality, which matches expectations set in the literature [15,27].

Conclusions
This work aimed to develop a methodology to generate machine learned material models while maximizing interpretability.To do this, we started with a known material model, the Gurson model for porous ductile metals, and systematically relaxed the associated assumptions.Relaxing the assumptions in this way made it so that newly generated GPSR models could correlate directly to the material's mechanics, as we controlled the assumptions in the training data.This systematic way of relaxing assumptions also allowed TL methods to be used by incorporating prior knowledge from the known model and information generated from previously generated models with fewer relaxed assumptions.The effectiveness of these methods was tested by relaxing the assumptions of void interaction (assumption 1) and void self-similarity (assumption 2) in the Gurson model.The methods were evaluated based on the convergence rate of the GPSR algorithm and the generated models, which were compared to the Gurson model to allow a direct correlation of additional terms to assumptions relaxed in the training data generation.
The convergence results on the assumption 1 dataset, figure 4, showed that the seeding algorithm had a significant impact on convergence as the inclusion of the Gurson model in initial GPSR populations reduced the fitness achieved in early generations.It was also shown that including the hyperbolic cosine operator was significant in achieving lower fitness values.For the fitness regularization algorithm, the penalty parameter, λ, was shown to perform better at a value of 5% than 20%, as it did not overly penalize models that violated the von Mises yield surface at VVF = 0 and was able to achieve similar levels of fitness compared to the other algorithms.The effect of the seeding algorithm on convergence was further shown on the assumption 2 dataset using two sets of seeds: models that abide by the von Mises constraint and models that do not.It was shown that the seeds that did abide by the constraint could achieve lower fitness levels on average, likely due to their lower complexity.
In producing physically reasonable models, the most successful algorithm was the fitness regularization algorithm, which penalized models that did not abide by the von Mises constraint.The final models were selected for analysis by down-selecting the models that abided by the von Mises constraint.Then, we showed multiple model selection methods, including best-fit and Pareto analysis, with advantages and disadvantages concerning fitness and interpretability.To balance these, we used GPSR as a feature engineering tool, selecting the most common terms and combining them into one model.Engineering models in this way resulted in two models, shown in equations ( 9) and (10) for assumptions 1 and 2, respectively.These models contained additional terms such as (σ vm /σ y )f and σ vm /σ y f that helped constrict the yield surface at larger values of von Mises stress and VVF.This additional constriction improved fit on training data over Gurson, especially in low stress triaxiality, which had previously been noted in many other studies [15,24,28].These additional terms can be attributed directly to the additional mechanics in the relaxed assumptions.
Comparison of the implicit fitness on the data when moving from the Gurson model, equations ( 1)-( 9), shown in figure 7, also showed improved performance in areas of low stress triaxiality and had an average relative difference of 1.42.The implicit fitness comparison moving from equations ( 9) and (10), shown in figure 10, resulted in a similar difference, with an average relative difference of 1.48.Assumptions 3-6, as listed in section 1.1, are still important since we have not considered damage evolution in the first two assumption relaxations.As further deformation is considered to the point of void coalescence, we can expect to see further differences in damage evolution.
The models generated in this study can be transferred to a wider variety of problems in ductile metal plasticity than the Gurson model due to fewer assumptions being present in their formulation.Due to the assumptions still being made in generation of these models, assumptions 3-6 in section 1.1, the accuracy of the model on other problems where these assumptions are not present will be limited.Modeling materials in this way allows a significant improvement in computational efficiency versus explicitly having to model the microstructure of a material as part of a larger structural model.For the computational cost you save, however, the computational cost for generation of GPSR training data does need to be considered.Implementation of these models in FE codes does require a state variable evolution equation for the VVF.In this case, equation ( 2) can be used as it still applies due to its derivation from conservation of mass.
The ability to impose TL techniques and interpret models in the ways shown here is due to our methods of systematic assumption relaxation with GPSR.While previous studies have used GPSR for material modeling applications, many lack a high level of interpretability due to a lack of understanding of the underlying assumptions in their experimental or numerical simulation data.Here, we control our assumptions explicitly to understand how they contribute to the final GPSR model forms, which maximizes the ability to interpret the models output from GPSR.Overall, we showed the effectiveness of this method in producing interpretable models and the potential of using TL techniques within GPSR.Future work will be dedicated to applying the techniques presented here to develop models under further relaxation of the remaining assumptions, and using these models to analyze the implications of the assumptions, ultimately deriving models that are both more interpretable and more widely applicable than their predecessors.

Figure 1 .
Figure 1.Process of systematically relaxing assumptions for generation of explainable material models.

Figure 2 .
Figure 2. (a) Data generated for self-similar void growth by FE modeling.(b) Individual load case taken from the self-similar void growth dataset with labels indicating time, t, and the time at the end of the simulation, t f .(c) Geometry of the RVE used to generate the datasets described in table 1 with applied triaxial boundary conditions via PBCs at the control nodes.Applied displacement ux+ controls the amount of applied displacement on the X+ face of the RVE.

Figure 3 .
Figure 3. (a) Acyclic graph representation and (b) operation list for the acyclic graph of the arbitrary model x0 + c0 sin(x1).Operations that do not contribute to the evaluation of the operation list are excluded.

Figure 4 .
Figure 4. (a) Convergence of the GPSR algorithm on the assumption 1 dataset for six tested configurations, where (b) is a zoomed-in view.Curves represent the average bestachieved fitness over ten trials.Fitness values of all algorithms are evaluated using equation (6).

Figure 5 .
Figure 5. (a) Pareto front for the 19 down-selected models.The plot in (b) represents a zoomed-in view to better observe differences in model fitness.

Figure 6 .
Figure 6.Engineered model for the assumption 1 dataset, equation (9), with complexity 19 and fitness 0.071 compared to the Gurson model.f is indicated on the colorbar.

Figure 7 .
Figure 7. Difference in implicit fitness calculated between the Gurson model, equation (1), and the engineered model, equation (9) on the assumption 1 dataset.Larger positive values illustrate better performance from equation (9).

Figure 8 .
Figure 8. Convergence of the GPSR algorithm on the assumption 2 dataset for two different seeding approaches.Regularized seeds are models from the assumption 1 Pareto front that abide by the von Mises constraint.Best-fit seeds are the best-fitting models from the assumption 1 Pareto front but do not have to abide by the constraint.

Figure 9 .
Figure 9. Engineered GPSR model for the assumption 2 dataset, equation (10), with complexity 26 and fitness 0.069 compared to the Gurson model.VVF is indicated on the colorbar.

Figure 10 .
Figure 10.Difference in implicit fitness calculated between equations (9) and (10) on the assumption 2 dataset.Larger positive values illustrate better performance from equation(10).

Table 1 .
Proposed models to address systematic assumption relaxations.

Table 2 .
GPSR hyperparameters used for training the relaxed assumption models.

Table 3 .
Algorithms tested on the assumption 1 dataset.