Reduced Order Numerical Modeling for Calibration of Complex Constitutive Models in Powder Pressing Simulations

Numerical simulations of different ceramic production phases often involve complex constitutive models, with difficult calibration process, relying on a large number of experiments. Methodological developments, proposed in present paper regarding this calibration problem can be outlined as follows: assessment of constitutive parameters is performed through inverse analysis procedure, centered on minimization of discrepancy function which quantifies the difference between measurable quantities and their computed counterpart. Resulting minimization problem is solved through genetic algorithms, while the computational burden is made consistent with constraints of routine industrial applications by exploiting Reduced Order Model (ROM) based on proper orthogonal decomposition. Throughout minimization, a gradual enrichment of designed ROM is used, by including additional simulations. Such strategy turned out to be beneficial when applied to models with a large number of parameters. Developed procedure seems to be effective when dealing with complex constitutive models, that can give rise to non-continuous discrepancy function due to the numerical instabilities. Proposed approach is tested and experimentally validated on the calibration of modified Drucker-Prager CAP model, frequently adopted for ceramic powder pressing simulations. Assessed values are compared with those obtained by traditional, timeconsuming tests, performed on pressed green bodies.


Introduction
Pressing and sintering are different manufacturing phases used to form parts starting from a powder.Pressing represents mechanical compaction of a powder in order to form the geometrically shaped solid body, called "green body", dimensionally close to the final part, which is obtained after a subsequent phase of sintering at high temperature.This production path is widely applied for manufacturing of various types of materials, including metallic materials through powder metallurgy [1][2], metal/non-metal composites [3][4], since offering larger flexibility when mixing components, ceramic materials [5], with high melting temperature, as sintering is generally performed at temperatures well below the melting one [6][7].Though it offers diverse advantages, such production is fairly sensitive to the selection of various parameters (e.g.pressing force, particle size distribution, sintering heating rate etc.).Even a small change in these parameters can affect properties of the final product, with some defects being evidenced only at the end of the production process.As an example, the lack of density homogeneity within the green body can cause different shrinking over the component during sintering, potentially violating the dimensional tolerances, or more severely can provoke the onset of internal cracks, resulting in component rejection.As a consequence, the production is usually relying on experienced operators to make final trial and-error adjustments [8].
Obvious practical advantages emerge in the above-outlined technologies by the employment of numerical methods for accurate simulations of diverse production phases.The development of these methods has been the focus of researchers, since the early works [9][10], considering only two-dimensional simulations, up to more recent studies dealing with multiscale three-dimensional models [11,15].The methodologies used to model the press and sinter powder manufacturing includes micro-mechanics [13][14][15], molecular dynamics [16] and continuum approaches [17][18][19].
Among the methodologies, continuum models have the benefit of shortest computing time, with the ability to predict attributes of interest, like density distribution, shape, residual stresses etc.
Constitutive modeling of material involved within pressing and sintering, required by continuum approaches is a challenging task.Many powders are formed as a mixture which can melt, react, diffuse and even alloy during sintering [8], clearly requiring multi-physics models with an elevated number of governing parameters.Similar complexity characterizes the constitutive description required for the pressing phase since the initial powder has fairly different behavior than the solid green body, formed during the process.During pressing green body is undergoing volumetric plastic deformation, changing also the elastic properties, so it is required to utilize elastoplastic coupling within plasticity models, with most governing parameters changing exponentially with relative density [19].
Calibration of such models is a difficult task, often involving a large number of experiments [20].As current praxis, the calibration of constitutive models used for powder pressing involves a series of destructive tests performed on green bodies [19,[21][22].An additional complication is encountered for certain powder mixtures where particles deform and undergo viscous flow to the die wall, effectively changing friction during the pressing.In such cases, a series of uniaxial compression tests are performed on a powder, with various levels of lubrication in order to obtain an estimate of friction coefficient value [8].
Clearly, such experiments are not suitable for the routine industrial applications.Furthermore, some parameters do not have clear physical meaning, so it is difficult to quantify them directly from the experiment [23].As a result, the accurate input data occasionally are not available, so approximations or simplified relations are used within the simulations.
An alternative advantageous strategy of the assessment of governing constitutive parameters is obtained by employing Inverse Analysis (IA) methodology.The present study aims at the development of a procedure for calibration of complex constitutive models used within pressing simulations, through IA methodology by means of data collected from pressing experiment only, without further experimentation on a green body.
The assessment procedure based on IA rests on minimization of a discrepancy function designed to quantify the difference between experimentally measured quantities and their computed counterpart.This minimization can be achieved by employing mathematical programming algorithms, like Trust Region Algorithm (TRA) [24], which implies calculation of first derivatives, numerically approximated by forwarding finite differences.In the present study, due to the complexity of considered constitutive model, the discrepancy function turned out to be rather convex, with unstable numerical simulations for certain parameter combinations.Such circumstance makes the TRA inefficient, since the value of the discrepancy function, or its derivatives, cannot be computed at every point in parameter space.
Alternatively, the minimization problem can be solved by means of soft computing techniques like Genetic Algorithms (GA) [25], capable of reaching a global minimum of no continuous function even with strong non-convexity.Such result is achieved through a sequence of function evaluations, typically larger by one order of magnitude than what is required by TRA.Significant time savings can be achieved if the model reduction is adopted for the non-linear test simulations needed to compute the value of discrepancy function.
In this paper, a novel procedure is proposed and investigated regarding the aboveoutlined problem of calibration.Developed procedure combines GA algorithm with test simulations performed by Reduced Order Model (ROM) based on Proper Orthogonal Decomposition (POD) and subsequent interpolation by Radial Basis Functions (RBF).This methodology has already been adopted in different research and described in details in [26][27].Here a modified version is developed and utilized, specifically addressing the laborious phase of "training", which is computationally the most demanding.This preliminary phase contains a fairly large number of simulations performed by the full numerical model (e.g. by FEM).In what follows a new strategy is proposed which utilizes ROM with controllable enrichment during the optimization by placing new training simulations within the zones with lower discrepancy function, where it is more important to have accurate calculations.Such strategy results in significant time savings as the recursive simulations within the minimization phase, are mostly performed by POD-RBF reduced model, gradually becoming more accurate throughout the minimization.
Developed procedure is experimentally verified with reference to the modified Drucker-Prager CAP model already employed by several authors for cold pressing simulations (see e.g.[19,[21][22]).With additional complexity adopted in this study, with respect to previous works, the number of parameters to identified equals 19.Resulting procedure turns out also to be numerically more stable as the fast POD-RBF model regularizes the convergence problems otherwise present with traditional FEM computations.The rest of the paper is organized as follows: Section 2 is devoted to a brief outline of the applied inverse analysis procedure for the assessment of material constitutive parameters in the present context; developed procedure which combines GA with specifically designed ROM is given in Section 3; within Section 4 validation of proposed strategy is presented in terms of comparative results considering values obtained by a novel procedure, and through more laborious experiments on green bodies; Section 5 contains conclusions regarding advantages and limitations of presented method, with potential future developments.

Experimental procedure 2.1. Constitutive model and the assessment of governing parameters through inverse analysis
Within present study a continuum approach is selected to model powder pressing phase, offering the advantage of modeling at one scale only.Modified Drucker-Prager CAP model was adopted as the constitutive model.Original Drucker-Prager (DP) yield criterion represents a pressure sensitive generalization of a popular von Mises criterion.DP yield criterion reads: where I 1 is the first invariant of the stress tensor, q is the equivalent von Mises stress, while d 0 and α are material parameters.Clearly, for α = 0, DP yield criterion reduces to von Mises.Such generalization is appropriate for materials with internal friction, like powders, as for these materials the pressure is influencing the onset of yielding.
To eliminate the deficiency of DP yield criterion, with material withstanding unlimited hydrostatic pressure, an additional elliptical yield surface is added (i.e."cap"), to form modified Drucker-Prager CAP yield criterion.Visualization of this double-surface yield criterion in the meridian plane is given in Fig. 1, with hydrostatic pressure reported as a positive value on the abscissa, typically adopted in soil mechanics, as opposite to solid mechanics convention.Additional CAP surface is defined by the following yield condition: Clearly, this addition introduces new constitutive parameters, here specifically: parameter R that controls the shape of the cap, parameter β used to define the transition between shear failure and cap surface, p a the evolution parameter related to the hydrostatic pressure yield stress p b .For a detailed description of modified form of Drucker-Prager yield condition, the reader is referred to [19].Further particularizations of the above model, adopted in the present study, are given as follows: • Relative density is related to the volumetric plastic strain, utilized as internal variable within the employed constitutive model.• Parameters are assumed to have an exponential dependency on relative density, namely for generic parameter P the relation reads: where 0 ρ and 1 ρ are fixed values of relative density, here specifically 0.5 and 0.9 respectively, while P 0 and P 1 are corresponding values of the parameter at these relative densities, and n is the exponent of this transition.Therefore, there are three values to be assessed related to each constitutive parameter: P 0 , P 1 and n. • The coefficient of Coulomb friction between the specimen and die wall is assumed as additional unknown and subjected to the identification.The evolution of parameter p b according to (3) is introduced in numerical model through hardening law, while for remaining parameters, namely: Young's modulus (E), Poisson's ratio (ν ), friction angle (α), cap eccentricity (R) and cohesion (d 0 ), additional implementation through ABAQUS user sub-routine [28] was required, in order to define these parameters as field dependent values (i.e.depending on volumetric plastic strain).
The above-outlined formulation represents a more realistic alternative to most of the previous studies, where diverse simplifications were used.For example, in [19] the influence of die wall friction was neglected, which in turn required a specific experiment with wall lubrication to minimize the frictional effect.In [29] cohesion and friction angle were considered as constants within compaction simulations, while in [30] it was assumed for Poisson's ratio.Clearly, more detailed description results in the more complex calibration process, which in the present study consists in the quantification of 19 parameters: 6 constitutive parameters as a function of relative density following relation is given by (3), and one friction coefficient, considered as a constant during the compaction.
Current praxis for calibration relies on the extensive experiments performed on green body apt to quantify the parameters, which are regarding the only single value of relative density (see e.g.[19,[21][22]29]).In order to quantify them over a wider range, it is required to perform a significant amount of experiments, with some values of relative density difficult to test.The advantageous alternative to this calibration is offered by the application of Inverse Analysis (IA) methodology, proposed in this study.
IA represents a synergic combination of experiments with numerical simulations and mathematical programming apt to provide a transition from experimentally measured quantities to required constitutive parameters.It is centered on the minimization of a "discrepancy function" that quantifies the difference between measured quantities and their computed counterpart.Such problem can be briefly formulated as follows.
Let u e be the vector of experimental data.The minimization with respect to the parameters in p, of a discrepancy function can be thus formulated if u(p) represents the measurable quantities related by test simulation (i.e."direct analysis") to the parameters sought as unknown variables: In the present study, digitalized force vs displacement curve, collected from pressing test, together with the value of radial stresses measured at two points along the height of compacting specimen are exploited to form a vector of measurable quantities.Details of used experimental configurations are given in Section 2.3.
The solution to the minimization problem defined by ( 4) directly provides values of sought constitutive parameters.Clearly, such formulation provides an important advantage as the material parameters are obtained straight from measured quantities within the powder pressing experiment, with no further experimentation needed on green bodies.In addition, considering the assumed relation (3), calibration through IA procedure results in the quantification of parameter value over a wide range of relative density.Obtaining numerical solution to the above minimization problem is therefore central within formulated inverse analysis.
The application of mathematical programming first-order Trust Region Algorithm, as the most efficient tool to find the numerical solution of minimization problems of the type given by ( 4), turned out to be inapplicable in the present case.Such circumstance is related to the large complexity of constitutive model subjected to the calibration, with unstable numerical simulations for certain parameter combinations.
To overcome this difficulty, specifically designed Genetic Algorithm (GA) was utilized in obtaining the numerical solution of above minimization problem (4).The minimization performed through GA is based on a sequence of "generations" over the parameter search domain.Main features of the algorithm designed and employed herein are outlined in what follows: • A set of N = 200 vectors of parameters is randomly generated, and for each of them, the discrepancy function is calculated.This set forms the initial "population", and the members are sorted in the ascending order of the corresponding value of the discrepancy function.
• A subsequent generation is formed by taking the first 2 members with the best value of discrepancy function and directly passing them to next generation as the "elite" members; the worst 50 members are subjected to "mutation", namely a random perturbation of their parameter vectors; the remaining members are subjected to a "crossover", where two members of this group are selected as "parents" and their parameter values are randomly combined in order to form two additional members with recombined parameters, to be passed to the next generation.• In a view of the above mentioned numerical instabilities, in the case of lack of convergency, the member is replaced by another one, with a slight modification of parameters so that the overall number of members in each generation is kept constant.
• The process is repeated until some of the specified stopping criteria are reached.From the last generation upon the termination of the optimization process, the best member is taken as the result of the minimization problem.The outlined scheme clearly is beneficial for considered minimization problem with respect to the TRA, as it is not requiring the discrepancy function to be, neither continuous nor convex over the parameter search domain.However, it involves a significantly larger number of simulations, which in a present context represents a sever burden and handicap if applied in routine industrial use.A remedy to this inconvenience, adopted here, consists in employing a Reduced Order Model (ROM) to perform test simulations.The outline of designed ROM which is incorporated within the optimization genetic algorithm is given in the following Section.

Reduced order model with controllable enrichment
In practical problems of parameter calibration through inverse analysis, a recurrent employment of the direct analysis (i.e.test simulations) is required.To make the resulting procedure more economical in terms of computing times, recourse is made to reduced order modeling, here specifically based on Proper Orthogonal Decomposition.Such a concept exploits the "correlation" between the computed test responses and is already applied in different contexts [27], [31][32].Details of the methodology can be found in [26], while in what follows an outline of the procedure is given, with particular emphasis on its incorporation within the optimization algorithm.
(a) Within the "search domain" in the space of parameters to be determined, N points are selected (called "grid nodes").(b) Each one of these N nodes is assumed as input of direct analysis simulating the test by utilizing full order model (i.e.FEM) and leading to vector (called "snapshot") collecting M test simulation quantities, here specifically displacements and corresponding forces required to form force-displacement curve for die compaction, and radial stress measurements.(c) Resulting N vectors are expected to be correlated since corresponding to the same mechanical system, with different input parameters, suggesting that a new "basis" can be computed in which axes with negligible components of vectors are dropped.(d) With the above model reduction procedure, each vector is now approximated by its "amplitude" vector in the new basis through a matrix generated by the eigenvalue computation of a symmetric matrix of order M. (e) Starting from any new parameter vector p, the corresponding vector a and the related test simulation quantities vector u, can now be computed with controllable accuracy and with much smaller computational effort by means of RBF interpolation among the previously computed responses with parameters by test simulations through full order model.
In the previous sequence, operations (a) -(c) are referred to as "training" of the model and can be time to consume.The selection of N grid nodes influences the accuracy of resulting ROM.The most controllable accuracy is achieved if each side of the search domain interval is subdivided into n intervals, but the number of such nodes grows exponentially with the number of sought parameters.This feature constrains the application of the outlined model reduction in the present context as the number of 19 sought parameters makes impossible the use or regular grid nodes distribution.To overcome this shortcoming, a scheme is proposed with employment of ROM with controllable enrichment, that is directly incorporated into the optimization genetic algorithm.The description of the proposed scheme in a form of sequence of operative steps is given here below.
a) A set of N = 200 parameter combinations is randomly generated over the search domain and for each of them, test simulation is performed in a force controlled regime by utilizing finite element model.From each analysis, a snapshot vector is formed collecting 100 displacements corresponding to 100 equidistantly spaced force levels from zero to maximum specified force.Measurements of radial stresses at two points along the specimen height are provided for the same 100 force levels.Account is taken through normalization for diverse orders of magnitude of different entries.5) -( 7) served to form reduced order model which can further be used to compute test response to an arbitrary set of parameters p through a matrix multiplication given by: where g is a vector of length N with entries computed through (6).e) N simulations performed within step (a) are further exploited to form the first generation for the optimization genetic algorithm.The members are sorted in ascending order of the discrepancy function and the second generation is formed by employing mechanisms summarized in Section 2.1.f) Subsequent generations are formed by performing all new computations through ROM, resulting in significant time savings as the computation given by ( 8) represents a straightforward matrix multiplication and it is therefore by several orders of magnitude faster than nonlinear FEM simulation.g) In the case of the replacement of the elite member due to the improvement of the discrepancy function, response to the test for that member has computed also through full order FE model for comparative purposes, to assess the error of reduced model.Should the error be larger than prescribed tolerance, additional FEM calculations are performed for best 50 members within a current generation with newly generated snapshots added as enrichment to those already computed within step (a).Refined ROM is now formed considering enlarged snapshot matrix by performing operations (b) -(d).
The outlined sequence is repeated until convergence criteria are met, specifically here regarding a number of generations during which no improvement in the discrepancy function was achieved.Such a scheme offers an important reduction in computing time as it is expected that most of the function evaluations are going to be performed by ROM.Further on, the gradual enrichments of ROM are performed within identified zones of the low level of discrepancy function, with model becoming more accurate particularly within the zone of interest.Proposed strategy is applied to the calibration of constitutive model outlined in Section 2.1 based on experiments performed on pressing of an alumina powder.Results are summarized in the following Section.

Adopted experiment and related numerical modeling
For validation purposes of the proposed calibration method pressing experiment of alumina powder is considered.The test is simulated by finite element modeling, utilizing the commercial code ABAQUS [28], adopting as material model Drucker-Prager CAP model with particularizations outlined in Section 2.1.
Developed calibration procedure is designed to exploit only data which can be collected from the pressing experiment.In order to make the more heterogeneous state of stress within compacting specimen, and therefore obtaining "rich" experimental data, in terms of being influenced by diverse constitutive parameters, two different configurations of die and punch are simultaneously utilized.Schematic representation of both configurations (further in the paper referred to as configuration 1 and configuration 2), together with the position of gauges used to measure radial stresses are visualized in Fig. 2. In a view of significant difference in rigidity between specimen and surrounding die, made of steel, both die and punch are modeled as rigid analytical surfaces.Between all surfaces a Coulomb friction unilateral contact is considered with coefficient treated as an unknown value, subjected to the identification procedure, together with other material parameters.To abbreviate computing time for test simulations, for numerical modeling of configuration 1, the axial-symmetry is exploited, with the model being two-dimensional.The lack of axial symmetry of configuration 2 enforced three-dimensional modeling, with half of the model being considered.The adopted finite element mesh for both models is visualized in Fig. 3.As specified in Section 2.2, 100 displacements and 100 radial stresses from both measured locations, corresponding to the prescribed level of pressing force are collected from experiments.The three curves, i.e. force-displacement, and two stress measurements, from each of the two experiments, served to form a so-called "residual vector", namely u e u(p), further used in discrepancy function within the minimization problem (4).Fig. 4 schematically shows how the entries of the residual vector are formed from considered curves.Additional press tests were performed utilizing the configuration 1 to produce cylindrical green bodies with different levels of relative density, specifically: 0.73, 0.78, 0.83, 0.88 and 0.93, to be further subjected to Brazilian test and Crush test.On the basis of these tests, through the procedure explained in [21], following constitutive parameters are calculated: cohesion (d 0 ), friction angle (α) and hydrostatic pressure yield stress (p b ).These values served as a reference to compare against those assessed through a novel procedure presented here.

Results and discussion
The solution of defined inverse problem on the basis of the procedure presented in what precedes, led to the parameter values listed in Table 1.  1. Simulated and experimentally measured curves are compared in Figs.5-7.6 shows that computed response, in terms of force-displacement curve, is reasonable close to the experimental one.On both figures also the curves with parameters used in the initialization are visualized, corroborating that the minimization process is done well even when starting from curves quite far from the "target" experimental ones.In the same manner Fig. 7 shows one of the radial stress measurements, that turned out to be representative of all other measurements, not shown here for brevity.
The above minimization result was achieved after 183 generations, each containing 200 members, hence involving 200×183 = 36600 function evaluations.The minimization was terminated by reaching stopping criterion of exceeding 40 generations with no improvement in the discrepancy function, with an overall reduction from its initial value of 24469.8 to a final one of 235.2.Within proposed procedure, the complete first generation and additional 19 enrichments, all including 50 simulations, were performed by full-order finite element model.A total number of time-consuming simulations, therefore, amounts 1150, which is less than 4% of the overall number of function evaluations.Considering that remaining function evaluations are computed on the basis of ROM through (8), such scheme contributes to the shortening of computing time by a factor of more than 30.It represents important time saving, as the pressing simulation for a three-dimensional model of configuration 2 takes about 10 minutes on a computer with i7 processor and 8GB of RAM.Reaching the same result through ROM is done in a real time.Selected parameters are compared against the values obtained through the procedure described in [21], as outlined in Section 2.3.The comparisons are visualized in Figs.8-10.
These results contribute to the conclusion that assumed relationship between parameters and relative density given by (3) seems to be general enough, resembling both exponential trend (for cohesion and hydrostatic compressive yield stress) as well as the linear one (friction angle).The accuracy of ROM is finally verified by performing additional comparison considering as inputs final set of parameters and performing computation both by full-order finite element model and by ROM.The comparison led to the results visualized in Fig. 11, evidence that there is the insignificant difference between the two curves.Such result corroborates the conjecture that the fast computations offered by reduced model are apparently not introducing any additional error within the identification procedure.

Conclusion
Material calibration procedure based on inverse analysis that is proposed and investigated in this study led to the conclusions which are briefly outlined in what follows.
Complex constitutive models with a large number of governing parameters are frequently employed in the simulations of diverse production phases of ceramic components.Calibration of such models through inverse analysis (IA) methodology turned out to be promising, offering larger flexibility with respect to conventional testing methods.With the example treated in this paper, it was demonstrated that IA provides an efficient transition from experimental data collected from pressing test to the constitutive parameters defined as a function of relative density.
The inverse analysis procedure inherently contains a numerical simulation of the test, which may be unstable when complex constitutive models are utilized.Such circumstance penalizes the use of efficient first-order (i.e.computing first derivatives only) mathematical programming algorithms for providing the solution to the resulting minimization problem within IA.Alternatively, genetic algorithms can be exploited, with increased computational effort, as these algorithms are involving a larger number of simulations.Reduced order modeling, as proposed herein, turns out to be crucial in order to make the computational time consistent with constrains related to the routine industrial applications.Proposed training scheme, different from previous studies (see e.g.[27,[32][33]) worked fairly effective in present context, by providing the enrichments of Reduced Order Model (ROM) throughout the optimization.The computational effort is therefore concentrated within the zone of parameters with a low value of the discrepancy function.Such strategy may be an advantageous alternative to the training of ROM according to "regular grid of nodes", for any situation where the number of sought parameters is high.This circumstance frequently occurs in complex multi-physics simulations utilized in ceramic productions.
Future promising developments current in progress concern the use of presented method to the calibration of constitutive models with even larger complexity used for simulations of other phases of production.The real challenging task, worth of investigation, is the application of such a scheme to multi-scale models.

11 ]=
in the previous step served to form snapshot matrix .According to POD theory (for details see[26]) by performing eigenvalue computation of matrix , a new basis is formed with orthogonal directions collected in a matrix Φ .The eigenvalues provide information exploited for truncation of the new basis, namely, by preserving only directions corresponding to the largest eigenvalues.In the present problem the drop by several orders of magnitude was evidenced for th eigenvalue with respect to the largest one, therefore truncated basis Φ )gathered first 10 directions.c) Matrix of "amplitudes", representing a projection of original snapshot matrix to a sub-space spanned by Φ ) RBF interpolation inverse multiquadric function was adopted which reads: coefficients, collected in matrix B, are computed through following matrix equation: considering as function arguments all parameter vectors generated in step (a).

Fig. 5 .
Fig. 5. Simulated and experimental force-displacement curve for configuration 1: computed with initial and with final parameter values

Fig. 6 .
Fig. 6.Simulated and experimental force-displacement curve for configuration 2: computed with initial and with final parameter values.

Fig. 9 .
Fig. 9. Comparison of friction angle (α) assessed through the novel procedure and on the basis of destructive tests performed on green bodies.

Fig. 10 .
Fig. 10.Hydrostatic pressure yield stress (p b ) assessed through novel procedure and tests on a green body.

Fig. 11 .
Fig. 11.The force-displacement curve for configuration 1 corresponding to a final set of parameters: comparison between FEM and ROM calculations.