Selecting the objective function during the inverse identification of the parameters of a material model of concrete

A BSTRACT . Selecting the correct objective function is the critical precondition for a successful optimization task. The validity of this condition is also required when optimization algorithms are needed for the inverse identification of the unknown parameters of nonlinear material models of concrete, where experimentally measured load-displacement curves can be conveniently applied. In such cases, the objective function expressions can be formulated as the difference between the functional values of the curves or via comparing the characteristic features, which comprise the area under the curve and also the maximum functional value. The proposed article brings a study of the influence of the different formulations of the objectives functions to achieving optimum in the inverse analysis using genetic algorithm. The numerical part of the study was performed in the ANSYS computational system with use of multiPlas library of elasto-plastic material models from which the model based on formulations of Menetrey and Willam was chosen. Sensitivity analysis.


INTRODUCTION
he application of advanced nonlinear constitutive models of building materials can be described as the approximation of mathematical modelling methods to the real behaviour of structures. Such efforts are, however, often complicated by the existence of a wide set of input parameters for such nonlinear models. Interest in the phenomenon of nonlinear behavior of concrete is in order of a wide range of use of this material in scope of many researchers. However, the construction of a correct constitutive relationship which is able to express this nonlinear behaviour for various types of loading appears to be problematic [1]. One of the basic problems which arise when formulating a material model for concrete is the different responses of the material to tensile and compressive load [2]. For this reason, several approaches are used for the mathematical description of the behaviour of concrete. One of these approaches involves the use of theory of plasticity [3]. Applications of theory of plasticity to the description of the behaviour of plain concrete can be found in the work of authors in [4], Willam and Warnke [5], Bazant [6], Dragon and Mroz [7], Schreyer [8], Chen and Buykozturk [9], Onate [10], Pramono and Willam [11], Etse and Willam [12], Menetrey and Willam [13], and Grassl [14]. The use of pure plasticity theory is not sufficient due to the gradual decrease in the stiffness of concrete due to the occurrence of cracks [1]. This problem can be removed when damage theory is used, i.e. by using an adequate damage model. However, as Grassl claims [15], independent damage models are not sufficient when the description of irreversible deformations and the inelastic volumetric expansion of concrete is required. Despite the above-mentioned limitations of both approaches, there are advantages to using both of them in mutual combination, and they can be combined further with other approaches formulated within the framework of nonlinear fracture mechanics. However, the use of combined material models results in a problem in real life in the form of the large amount of parameters which need to be known for the selected material model before the launch of the numerical simulation itself. Unfortunately, not all data regarding these parameters, which can be both mechanico-physical and fracture-mechanical in nature, may be available in advance. This problems can be resolved with use of other advanced mathematical approaches like optimization analysis. The varied spectrum of possibilities characterizing the use of optimization methods includes, among other options, the above mentioned inverse identification of the unknown parameters of the nonlinear material models utilized in numerical analyses performed with the finite element method. The optimization algorithms exploited in the inverse analysis of unknown material parameters, described within references [16,17], constitute a counterpart to methods based on the training of artificial neural networks as discussed by Novak and Lehky [18]. However, both in cases where optimization is applied to the identification problem and during any classic use of the optimization methods presented in [19], the decisive factor to support a successful optimization process consists in selecting the appropriate algorithm and correctly formulating the relevant objective function. The actual need of such a function becomes even more prominent in identification using optimization modules implemented within ANSYS Workbench [20], where the definition and computation of the objective function value have to be performed with an external program or script. The task embodying the inverse identification of unknown material parameters consists in utilizing the experimentally measured curves that characterize the relationship between the load L and the deformation d (L-d curves). During the actual identification, this reference pattern is compared with the L-d curves produced by the nonlinear numerical simulations within the corresponding experiment. The basis of the comparison then rests in calculating the similarity ratio, which is represented by one or more numerical values and also prescribes the objective function. With respect to the formulation of the objective function, the optimization task is, in a given case, defined as the minimization of the similarity ratio. For the discussed purpose, it appears advantageous to employ the RMSE (Root-Mean-Square Error) ratio, an instrument that, according to [21], enables us to compare the differences between the values measured and those generated via a mathematical model; in this context, the authors of reference [22] analyze the application of the RMSE ratio within disciplines such as meteorology, economics, and demography. Considering the shape of the L-d curves, it is then possible, as shown in study [23], to exploit them in comparing the value of the surface below the loading curve with the maximum load value. Importantly, if the second one of the described variants is used, we also have to select a correct and robust algorithm to facilitate the optimization including multiple objective functions; this problem can be further encountered in the computation of more RMSE ratios for partial sections of the curve, whose positive impact is embodied in the analysis of the individual parameters' sensitivity to specific sections of the loading curve. With respect to the above-outlined conditions, the present article examines the effect exerted by different formulations of the objective function in a given identification task. For the identification proper, we chose the L-d curve measured during a three point bending test on a notched concrete beam, according to [24]. The numerical simulation of the experiment was performed with ANSYS via a nonlinear, multi-parametric material model of concrete adopted from the multiPlas library [25]. Generally, the paper aims to describe the applicability of the above-mentioned possibilities of formulating the objective function; the computation of the individual options was enabled by scripts created in Python.

INPUT DATA
n order to analyze the suitability of the selected objective functions, we chose one L-d curve associated with the set of fracture tests published by Zimmermann et al. [24]. The specimen, a notched concrete beam manufactured from class C25/30 concrete and having the length l equal to 360 mm, height h of 120 mm, width w corresponding to 58 mm, and notch height of 40 mm, was configured for three point bending test; the vertical deformation d was measured in the middle of the span of 300 mm at the bottom side of the specimen. The cited article presented experimental and numerical research where the identification procedure based on utilization of neural network was used. The identification I process was aimed to finding values of basic mechanic-physical and fracture-mechanical properties of concrete specimens: modulus of elasticity Ec, tensile strength ft and fracture energy Gft. The resultant statistical values of these properties were: modulus of elasticity E c = 38,9 GPa, tensile strength f t = 2,76 MPa and fracture energy G ft = 215,6 J/m 2 .

GEOMETRY AND MESH OF THE COMPUTATIONAL MODEL
he computational model and the nonlinear computation of the fracture experiment task were performed using ANSYS Workbench. To reduce the computing time, the geometry of the computational model was covered with a mesh of 2D finite elements (PLANE182) having the edge length of 6 mm, and the problem was solved as a plane stress task with the element thickness of w = 58 mm. The real support provided by steel bearings was, in the computational model, idealized via strain boundary conditions, where the vertical deformation was prevented at the location of the support, and the horizontal deformation was -with respect to the solvability of the task -restrained in the middle of the span at the upper side; such an arrangement then corresponded to the position of the introduced load. The notch was modeled in a simple manner, using a pair of parallel lines having a common node at the top of the notch.

MATERIAL MODEL
o ensure both the nonlinear behavior and the identification of the corresponding material parameters, we selected a material model from the multiPlas database; this model is based on the plasticity surface derived by Willam and Warnke [5] and then modified by Menétrey [26], Menétrey and Warnke [13], and Klisinsky [27]. Further, the applied material model represents the evolutionary branch of nonlinear material models of concrete that exploit plasticity theory combined with instruments of nonlinear fracture mechanics to form a single model. This product then ranks among the group of models with non-associated plastic flow rule and is formulated such that it considers the invariants of the stress tensors and the deviatoric stress tensors; thus, the plasticity surface edges are softened, and the definition fidelity of the nonlinear behavior of concrete improves [10]. The formula for the yield surface as found in the programme manual [10] has the following form: with the elliptic function r(θ,e) developed by Klisinsky [27] on the basis of Willam and Warnke´s findings [5] and where Ω(σ,κ) is a hardening/softening function with a work-hardening law and A, B, C, D are model parameters containing basic mechanico-physical properties of concrete (uniaxial compressive strength fc, uniaxial tensile strength ft and biaxial compressive strength f b ) whose form is given by the following relationships: The Eq. (1) also contains coordinates of Haigh-Westergaard cylindrical space, where χ represents the height, ρ the radius and θ the azimuth. These coordinates are functions of the above mentioned invariants of stress tensor and can be written in the following manner: The formula for the plastic potential takes the following form: T where parameters X and Y are again related to uniaxial compressive strength fc, uniaxial tensile strength ft, biaxial compressive strength fb and the so-called dilatancy angle ψ. and where As regards using the finite element method, the model utilizes the concept of smeared cracks [28], and, thanks to exploiting Bazant's crack band theory [29], it does not exhibit a negative dependence on the size of the finite element mesh. With respect to the above-described material model formulation consisting in a combination of several theoretical approaches, the total of 12 mechanico-physical and fracture-mechanical parameters were required to ensure the functionality of the model; a description of these parameters is presented in Tab. 1 [30].

INVERSE IDENTIFICATION
he lack of knowledge of the input parameter values constitutes, together with the continuous theoretical development in the given field, a basic problem affecting the actual performance of advanced nonlinear material simulations of concrete structures. However, such knowledge deficiency can be advantageously eliminated via inverse analysis and experimental research. Yet the use of optimization techniques to identify the discussed parameters T still poses questions concerning the choice of a correct formulation for the relevant objective function. The inverse identification within this paper was carried out with an optimization module implemented in ANSYS Workbench; using such a computing system nevertheless required us to create external scripts in Python, and these were called during the batch calculation to compute the relevant objective function values. The whole process of inverse identification consisted in sensitivity analysis of the input material parameters to shape of the L-d curve and optimization itself which was used for minimization of the difference between the numerical and experimental L-d curves. The optimization algorithm was used for varying values of the material parameters and the parameters belonging to the curve with the lowest value of the objective function could be then considered as the sought material parameters. The calculation was performed automatically via ADPL macro that prepared geometry and mesh of the computational model, set up the material model with appropriate values of the material parameters, solved the task and called external Python script for calculation the of the objective function.

Description of the Selected Objective Functions
The basic objective function giving the difference between two curves was embodied in the RMSE ratio. The calculation of this function could not be performed directly, because the distribution of points on the reference and numerical curves was invariably different due to the varied runs of the solver. The mapping of the points on the numerical curve according to the reference curve was, within the script, resolved via linear interpolation. After aligning the points on the curves, we calculated the RMSE ratio according to the formula  

Sensitivity Analysis
The actual identification of the material model parameters was invariably preceded by a sensitivity analysis aimed at mapping the space of the design variables and determining the sensitivity of the individual material parameters to the value of the objective function. For each identification task, we conducted 250 simulations, and uniform covering of the design space was ensured via the LHS method. The sensitivity of the individual material parameters to the output parameters was expressed using the Spearman correlation coefficient r s . The sensitivity analysis for the option with one RMSE optimization function showed that the highest sensitivity rate could be found in the elasticity modulus E, uniaxial tensile strength ft, and specific tensile fracture energy Gft. The high sensitivity rate of these parameters can be explained by the very basis of the examined problem: the simulated task is one with tensile bend. Very interesting results were obtained from the sensitivity analysis involving the subdivision of the material parameters into L-d curve sectors represented by 5 values of the RMSE ratio. In the given case, we identified that the first sector is utterly dominated by the sensitivity to the elasticity modulus E; the second portion is governed by the tensile strength ft; and the last part exhibits major influence of the specific tensile fracture energy Gft.
In the analysis of the parameter sensitivity to the values ΔA Ld and ΔL max , only the sensitivities to the tensile strength f t and the specific fracture energy G ft were revealed, which nevertheless cannot be considered correct with respect to the character of the task.

Optimization
The actual identification of the material parameter values was carried out via direct optimization using a genetic algorithm, and, considering the results of the sensitivity analysis, it was performed in a space of three variables. The reduced design vector then assumed the form As already indicated above, the identification in its entirety was performed three times; in the second and third tasks, however, two optimization variants were carried out. The actual process comprised only a minor formulation change, where the first phase involved minimizing the objective function values, and the second one consisted in seeking the zero value of the relevant objective function.

RESULTS
o present and compare the results, we selected the RMSE ratio, which was then computed in all calculation options. Within the initial identification task, the optimization algorithm generated 76 design vectors in total, and the minimum exhibiting the value of RMSE = 144.92 N was achieved with the 33 rd iteration. The second optimization task, or, more concretely, its variant that sought the minima of the RMSE ratio values within the L-d curve sectors, eventually generated 58 design vectors, and the minimum having the total value of RMSE = 160.39 N was obtained during the 10 th iteration. The second version of this identification task produced the minimum with the total value of RMSE = 136.92 N in the 165 th out of 188 iterations The material parameter identification conducted with the objective functions defined as the differences between ΔALd and ΔL max proved the applicability of the given manner of formulating the objective function, though only at the cost of the longest computational time (compared to the related options). In the minimization variant, we obtained the minimum exhibiting the total value of RMSE = 143.13 N during the 85 th out of 122 iterations, while the associated version provided the minimum at RMSE = 178.39 N in the 254 th out of 415 iterations. Comparative diagrams representing the resulting L-d curves and the history of the objective function values within the optimization process are shown in Figs. 1(a) to 1(c).

CONCLUSION
he results obtained within the presented research indicate that a successful identification can be performed with both an objective function defined as the RMSE ratio and an objective function formulated as the differences between the characteristic signs of the L-d curves: the surface under the ΔALd and ΔLmax curves. However, it can be also claimed that the best optical agreement was achieved with the objective function defined as the RMSE ratio along the entire L-d curves. By extension, the outcome of the research then points to the fact that the actual formulation of the objective function is influenced by whether or not the given function is only minimized or a zero value is sought. Even though the use of the variant seeking the zero value of the objective function was accompanied by increased computational time requirements, both of the above-characterized methods can be recommended for practical computation; now, however, it is also necessary to consider the fact that, besides the selection of the correct objective function, the choice of a suitable algorithm constitutes a major, decisive aspect within the discussed procedures.

ACKNOWLEDGEMENT
he research presented within this paper was supported from project GA14-25320S titled "Aspects of the Use of Complex Nonlinear Material Models" and guaranteed by Czech Science Foundation. The authors also wish to acknowledge the assistance provided by the project No. FAST-J-16-3562, "Implementation of Material Models of Concrete in ANSYS and their experimental verification", associated with the specific university research programme of Brno University of Technology. T T