Identification of Concrete Material Model Parameters Using Optimisation Algorithms

The application of nonlinear material models of concrete within numerical simulations focused on the design of safer and more economical protective concrete structures is currently the subject of investigation of many scientific researchers. However, one basic problem related to the nonlinear modelling of concrete is that very often there is a lack of knowledge about the material model parameters whose values must be defined. The solution to this problem can be in what is termed as inverse parameter identification, an approach which is presented in this paper. Specifically, the material parameters of the Continuous Surface Cap Model for concrete are identified within this paper using optimisation algorithms. The subsequent comparison of parameter identification results with experimental data shows the efficiency of the presented approach.


Introduction
Complex nonlinear material models of concrete implemented in sophisticated computing systems based primarily on the implicit or explicit finite element method [1][2][3] are currently the most powerful tools enabling the modelling of the real nonlinear behaviour of concrete within static and dynamic numerical simulations [4][5][6][7][8][9][10][11]. Such simulations can be focused on the design of safer and more economical protective and military concrete structures. However, the application of nonlinear material models of concrete within continuum mechanics tasks causes some problems which might arise and with which the designer must somehow cope. The most fundamental and largest

Modelling of Concrete
Within this paper, concrete was modelled via a nonlinear material model known as the Continuous Surface Cap Model [17,18]. This model can be found in the library of material models implemented in LS-Dyna software.
The Continuous Surface Cap Model is based on a yield surface which is defined as a function of three stress invariants according to the following equation [19,20]: where I1 is the first invariant of the stress tensor, J2 and J3 are the second and third invariants of the deviatoric stress tensor, ℜ(J3) is the Rubin strength reduction factor and κ is the cap hardening parameter. The yield surface is composed of two parts, these being the hardening compaction function Fc(I1,κ) and the shear failure function Ff (I1). The mathematical expression of the hardening compaction function is given by the following equations: Identification of Concrete Material Model Parameters Using Optimisation Algorithms 35 with the equations: where R is the cap aspect ratio. The shear failure function is defined by the equation: where α, β, λ, and θ are material constants usually determined on the basis of triaxial compression tests. Within the yield surface, the shear failure function and hardening compaction function are combined using a multiplicative formulation which allows their combination to be continuous and smooth at their intersection. Within its formulation, the model allows the effect of strain rate on the resulting stress state to be taken into account. However, this capability of the model can be neglected within the computations by appropriate setting of the model parameter IRATE (IRATE = 0: Rate effects turned off; IRATE = 1: Rate effects turned on). If the parameter IRATE is equal to zero, the response of the model is static, independent of loading rate. It follows that the Continuous Surface Cap Model can be used in dynamic (but also quasi-static or static) numerical simulations. The parameter IRATE equal to zero has been used within this study / research, so the response of the model was always static. It is also important to mention that the model includes an algorithm for limiting the dependence of results on the size of the finite element mesh.
The Continuous Surface Cap Model is implemented in LS-Dyna software in two versions, specifically the *MAT_CSCM version and the *MAT_CSCM_CONCRETE version [16]. The general *MAT_CSCM version of the material model was used within this study / research because this model version contains a lot of parameters and is therefore suitable for identification purposes. A total of 25 parameters (material constants) of the *MAT_CSCM version were identified in order to obtain the most realistic response from the model, i.e. a response that is as close as possible to that described by experimental data. Descriptions of the identified parameters are provided in Tab. 1 [16] along with the units used. The material model version *MAT_CSCM_CONCRETE was not used in the identification task performed within this study / research because this version requires the definition of only 3 parameters (mass density, the uniaxial compressive strength and the maximum aggregate size), based on which the other parameters are automatically generated. As a result, the postpeak response in particular of this model version can only be influenced very slightly. However, the results of the *MAT_CSCM_CONCRETE version were used for comparative purposes.

Experimental Data
For the parameter identification performed, use was made of experimental data obtained from the evaluation of direct tensile test results presented in [21]. Direct tensile test is one of possible tests which are used for investigation of tensile behaviour of plain concrete. The experimental data were represented by a force-extension curve which characterised the nonlinear behaviour of plain concrete test specimens during direct tensile loading (see Fig. 1).
According to [21], the dimensions of the plain concrete test specimens used during the tests were 305 mm × 60 mm × 19 mm (length × width × depth for the critical cross-sectional area). The experimental data were measured from an 85 mm gauge length of each test specimen. The 28 days uniaxial compressive strength of the used concrete was 44 MPa. The size of aggregate varied, however, the maximum aggregate size was 10 mm. During the direct tensile loading, concrete test specimens were stretched at a constant loading velocity. Loading was slow, so the response of test specimens was static. Fig. 1 shows the typical behaviour of plain concrete in direct tension. The graph depicts the linearly elastic behaviour of the concrete before it reaches the ultimate uniaxial tensile strength with subsequent nonlinear tensile softening.

Identification of Concrete Material Model Parameters
Using Optimisation Algorithms 37

Nonlinear Numerical Simulations
Within this paper, nonlinear numerical simulations were performed in the LS-Dyna software based on the explicit finite element method. Settings for the calculation, computational model and parameter values of the material model were always entered into the LS-Dyna input file which, together with experimental data, formed the input data for the optiSLang programme used for parameter identification.

Fig. 2 The finite element model of the measured region of the test specimen
The computational model used for this test was simplified in comparison with real direct tensile tests. Experimental data were only measured on the aforementioned 85 mm long region of the test specimen. This means that the extension in Fig. 1 corresponds to the extension of that 85 mm long region and the tensile force corresponds to the tensile force needed to obtain the given extension of that region. This is why it was the only part of the test specimen that was modelled. Explicit 3-D structural finite elements were used for the modelling. The geometric model was discretised by a regular mesh of finite elements. The size of the finite elements has been chosen so that the length of the time step in explicit algorithm was not too small. This had resulted in a reasonable calculation time and, therefore, the task could be solved as a whole without the use of symmetry. In terms of boundary conditions, linearly increased displacements over time in the X-axis direction were prescribed for the nodes of both finite element model bases. These displacements were applied so that tensile loading could be imposed at a constant velocity. Explicit finite element algorithm allowed to solve the given task without supports while maintaining computational stability. This was taken into account and no supports were applied to the finite element model. This enabled the proper transverse shortening of the model due to Poisson's ratio. The finite element model of the measured region of the test specimen is shown in Fig. 2.

Parameter Identification Process
The whole process of parameter identification was performed via the optiSLang programme and consisted of three parts: 1. Sensitivity analysis 2. Global optimisation 3. Local optimisation

Sensitivity Analysis
Within this first part of the whole parameter identification process, the sensitivity of identified material parameters to the defined reference response was analysed. The reference response was represented by individual points which defined the shape of the used force-extension curve. The main aim of the sensitivity analysis [22,23] was to determine the minimum scope of the design vector which contained the identified material parameters. Another aim was to modify the range of variability for the individual material parameters contained in the design vector. Boundary values of the initial range of variability for each identified parameter were obtained on the basis of test calculations. The Advanced Latin Hypercube Sampling (ALHS) method [24] was used to perform the sensitivity analysis. A total of 250 random realisations of the design vector were generated through this method. The number of random realisations covered the given design space sufficiently. However, the first realisation corresponded to randomly user-selected parameter values, for which the force-extension curve is illustrated in Fig. 3 together with the experimental data for comparison. The parameter values used for the first realisation respected the initial ranges of variability for individual identified parameters.
The results obtained via the ALHS method showed that only 9 out of a total of 25 identified material parameters significantly affected the resultant shape of the numerically-simulated force-extension curve. The influence of parameters on the forceextension curve was indicated by the Relative Frequencies of Importance (RFI), which were nonzero for significant material parameters (see Tab. 2). For the subsequent global optimisation, the original design vector which contained all identified material parameters given in Tab. 1 was reduced to a design vector, which only contained the aforementioned 9 parameters that significantly affected the resultant shape of the numerically-simulated force-extension curve. Specifically, the reduced design vector had the following form: The values of the pre-optimised material parameters and objective function (OBJ_FUNC) obtained from the best random realisation generated by the ALHS method are presented in Tab. 2.

Identification of Concrete Material Model Parameters
Using Optimisation Algorithms 39

Global Optimisation
Within this second part of the whole parameter identification process, optimised material parameter values were sought that would enable the numerical simulation result to approximate the experimental data very well. In other words, optimised material parameter values were sought that would minimise the value of the objective function. The global optimisation process was therefore based on minimising the objective function [15]. The objective function used within this paper can be formulated by the following equation: where ysim,i are the force values obtained from the appropriate numerically-simulated force-extension curve at the given displacements and yref,i are the force values obtained from the used experimentally-measured force-extension curve at the same displacements. As already mentioned, global optimisation only involved those material parameters which were included in the reduced design vector. Other parameters were defined by constant values (deterministically) from their modified ranges of variability. For the purposes of global parameter optimisation, an optimisation algorithm generally known as an Evolutionary Algorithm (EA) [15] was used. An EA is one of the optimisation approaches which exploit processes inspired by biological evolution, including, for example, reproduction, mutation and recombination. The five best random ALHS method realisations were used as a starting point for the EA. The history of objective function values during the run of EA is showed in Fig. 4. The curve in Fig. 4 describes the evolution of objective function values depending on the EA generations (optimisation designs). The best EA generation was the generation number 112 in which the objective function got the minimum value and the material parameters got the optimum values. It can be seen in Fig. 4 that the optimisation calculation continued even after the best EA generation was reached. However, once the algorithm evaluated that the results are no longer improving, it automatically terminated the optimisation calculation.
The optimised material parameter values from the best EA generation are presented in Tab. 2, including the minimised objective function value.

Local Optimisation
Within this third part of the whole parameter identification process, just as in the case of the global optimisation, optimised material parameter values were sought that would minimise the objective function value. Of course, the objective function used within the local optimisation corresponded to the objective function described in Eq. (9). The aim of the local optimisation was to search the vicinity of the global minimum for the purpose of its refinement. Just as in the case of the global optimisation, only the material parameters included in the reduced design vector were optimised while the other parameters were defined deterministically. The local parameter optimisation was performed using the direct optimisation algorithm known as the Simplex method [15]. The Simplex method belongs to those iterative algorithms that are carried out systematically to determine the optimal solution from a set of feasible solutions. The best EA generation was used as the starting point for the calculations performed using the Simplex method. The history of objective function values during the run of Simplex method is showed in Fig. 5 where it can be seen that the best Simplex method generation was the generation number 202. In this Simplex method generation, the objective function got the minimum value and the material parameters got the optimum values. As in the case of the EA, optimisation calculation continued even after the best Simplex method generation was reached but once the method algorithm evaluated that the results are no longer improving, it automatically terminated the optimisation calculation.
The optimised material parameter values from the best Simplex method generation are presented in Tab. 2, including the minimised objective function value. It can be seen from Tab. 2 that the Simplex method provided the best results because the objective function value obtained via this method is lower than the objective function values gained for the Evolutionary Algorithm and ALHS method.  Fig. 6 shows a comparison of the experimentally-measured force-extension curve with the force-extension curve obtained from the numerical simulation in which the resulting optimised parameter values of the Continuous Surface Cap Model from local optimisation (Simplex method -the best results) were applied. It can be concluded from Fig. 6 that the material parameters of the Continuous Surface Cap Model were optimised very accurately because the numerical simulation result ensures very good approximation of the experimental data.

Results
For comparison, the numerically-simulated force-extension curve for the material model version *MAT_CSCM_CONCRETE is illustrated also in Fig. 6. The three necessary parameters of this model version were defined by the values 2.4 × 10 −9 Mg/mm 3 for the mass density, 44 MPa for the uniaxial compressive strength and 10 mm for the maximum aggregate size according to the experimental data. It can be seen in Fig. 6 that the numerically-simulated data for the model version *MAT_CSCM_CONCRETE exhibit significant differences compared to the experimental data in the region of post-peak behaviour. Consequently, it can be concluded that the *MAT_CSCM version is more suitable for obtaining the best approximation of the experimental data than the *MAT_CSCM_CONCRETE version. However, in terms of the quantity of the material parameters, the *MAT_CSCM model version usually requires the use of the parameter identification process presented in this paper.
For another comparison, the numerically-simulated force-extension curves for optimised material parameter values are illustrated in Fig. 7. The curves were obtained for different loading rates and, for all curves, the material model parameter IRATE was equal to zero (rate effects turned off). The solid curve in Fig. 7 corresponds to the loading velocity of 0.537 mm/s that was used for all calculations performed in this paper. The dotted curve corresponds to ten times faster loading than is the case with the solid curve and the dashed curve corresponds to hundred times slower loading than it is in the case with the solid curve. It can be seen in Fig. 7 that for material model parameter IRATE = 0, the curves for all tested loading rates are practically identical

Identification of Concrete Material Model Parameters
Using Optimisation Algorithms 43 with very small differences. This confirms that if the parameter IRATE of the model is equal to zero, the response of the model is always static, independent of loading rate.

Conclusion
This paper was focused on performing the identification of the material parameters of the Continuous Surface Cap Model using optimisation algorithms. This was carried out utilising an experimental force-extension curve obtained on the basis of evaluating the results of direct tensile tests performed on specific plain concrete test specimens. The obtained results from the parameter identification process showed that the Continuous Surface Cap Model is a suitable tool to describe the nonlinear behaviour of real plain concrete in direct tension. They also indicated that, in cases where the parameters of the material model are conveniently selected and entered, we can obtain a very good agreement between the experimental data and numerical simulation. Importantly, this claim was proved by the outcome of the numerical simulation performed using the resulting optimised parameter values of the material model from the Simplex method. The discussed results approximated the used experimental data very accurately. This paper has also shown that with the help of powerful optimisation tools it is possible to obtain very good fit for the parameters of a nonlinear material model intended for concrete modelling. Advantageously, the product of the parameter identification process performed can be exploited for further research concerning the nonlinear numerical response of protective and military concrete structures and, of course, general concrete structures also.