Identification of the hyperelastic constitutive model constants of a rubber-like material

In this article, we present a possible strategy for solving an inverse problem with modeFRONTIER.
The solution of an inverse problem is often required in engineering situations, more often than one would expect. Most of the times, the values of quantities which influence the behavior of the system under exam are unknown and cannot be determined easily or directly.
Let us imagine, for example, a situation where we suspect an embeddedcrack in a large concrete dam. Obviously, only non destructive tests can be performed in this case in order to check if the position and the dimensions of the crack may eventually compromise the structural efficiency.

Fig. 1 - The bell-shaped vibration absorber considered in this example; an isometric view (a) and a vertical section (b).

In this situation dynamical tests can be performed on the dam (accelerometers are usually employed to monitor the response of the structure subjected to known sources of vibrations); a virtual model can be subsequently built and tuned using the collected data and used to predict the crack position and dimensions.
Generally, whenever some observed data of the system response are available, we can try to solve an inverse problem with the aim to find appropriate values for the unknown parameters.

The expression “inverse problem” comes from the fact that these problems are formulated in an “inverse way” with respect to the traditional “direct” approach to engineering problems, which usually require to determine unknown effects starting from known causes.
The efficient solution of an inverse problem can be a very difficult task, mainly because it is necessary to reproduce the system behavior in an accurate way, to have an exhaustive dataset of observed data and, finally, to have algorithms able to drive the search towards reliable solutions.
Sometimes, the solution of an inverse problem is not unique; many different sets of the input parameters can lead to similar outputs. For this reason, the solving procedure has to be as much robust as possible in order to capture all these possible solutions.

In this article we consider a simple example, where the objective is to identify the constitutive law parameters for a rubber-like material starting from some experimental observations. A bell-shaped vibration absorber device is considered (see Fig. 1): it is made of two steel components, the support base and the bell, connected by a rubber ring. This kind of devices is often used for the connection of vibrating machines (engines, rotating unbalanced mechanisms, working machinery, etc…) with the underlying structures, in the attempt to reduce the transmitted vibrations and absorb energy, thanks to the viscoelastic properties of the rubber.
In this example, the goal is to characterize the rubber-like material by means of a static non-destructive test. The experimental investigation has been performed directly on the dumper, increasing the vertical load, applied in a static manner, and measuring the corresponding vertical displacement of the bell.

The finite element model
The ANSYS Workbench 10 software has been used to build an axisymmetric model of the dumper.
The steel support is not rigorously axisymmetric, but this simplification is considered here to be licit, being the vertical deformation of the bell and the rubber ring not significantly influenced by the support shape. Moreover, when compared to an analogous three dimensional model, an axisymmetric model has the desirable advantage to require much less computational effort to be solved.
A three term Mooney-Rivlin model is adopted to describe the rubber hyperelastic behavior; other choices are of course possible. The strain energy function can be written as follows:
 
where I1 and I2 are the deviatoric strain invariants and J is the ratio between the actual and the initial volumes. C10, C01, C11 and D are the free parameters which completely characterize the material response (see Fig. 2) of the rubber.

Inverse problem definition
The objective is now to find appropriate values for C10, C01, C11 and D (hereafter referred to as the input parameters) such that the difference between the virtual and the real response is minimized.

Typical stress-strain curves obtained with different stress conditions for a rubber-like material

This last objective has been summarized in the following equation:

  where ym and yv represent the measured and the virtual vertical displacement of the bell, being N the number of available observations.
The choice of using just one target, as reported in the above equation, is mainly motivated by two considerations. The first one, is that the problem is reduced to a mono-objective optimization instead of a multi objectives one, with obvious simplifications; the second reason is that, from an engineering point of view, it is more interesting to have the best solution “in average”, rather than a solution that tries to exactly match all the experimental data, which could be affected by measurement errors.
It is clear that the usual data fitting techniques [3] used to find the strain energy parameters starting from experimental data, cannot be used in this case. These techniques actually require stress-strain measures obtained in very simple conditions (typically uniaxial, biaxial, pure shear and volumetric test conditions) which are not available in this context. For this reason the use of an optimizer, able to drive a simulation software is mandatory to identify the strain energy constants.
It is interesting to note that also other mechanical properties of the rubber material could be identified with similar procedures. The viscoelastic behavior for example, which could be added and coupled with the hyperelastic response, is usually modeled by a so-called Prony series whose terms could be considered as additional input variables.
The solution of this new inverse problem obviously would require other experimental tests, able to detect the dynamic response of the dumper.

Problem set up in modeFRONTIER
The first step consists in the construction of an appropriate work-flow, which has to include four input variable nodes (C10, C01, C11 and D), six output variables, one for each calculated vertical displacement of the steel bell at different load steps (S1 up to S6), the ANSYS Workbench node, the DOE and scheduler nodes, the logic end and the objective node containing the target function described above. A constraint which limits the initial shear modulus to positive values has been added.

It is also useful to introduce a vector input variable node to store the constant experimental data (see Fig. 3). The NSGA-II optimization algorithm has been chosen, considering that all input variables are continuous and that robustness is a desirable feature in this case. Once a solution has been found it could be subsequently refined with a gradient based algorithm, such as the Levenberg-Marquardt, in view of the peculiar form of the target function.

Fig. 3 - Deformed shape corresponding to the best solution found	Fig. 4 - The objective function value plotted versus the design ID. The chart also shows the moving average of the last 30 designs (the green line). It can be seen that NSGA-II tends in average to improve the solution through the generations.

The DOE consists in 30 design randomly generated. The initial population has been expressly generated with a relative large number of design being height the expected number of error designs.
Run the parameter search and interpreting results
Once the work-flow has been completed the parameter search can run. As mentioned above, a large number of error designs is present.

The main reasons of this fact are a non appropriate choice of the material parameters and a non optimal set up of the solver, which sometimes is not able to overcome convergence difficulties.
However, it is interesting to note that the error designs tend to decrease with NSGA-II generations; this is a very notable feature of genetic algorithms, which are able to automatically drive the solution search towards feasible regions.

Table 1: The best solution found using NSGA-II algorithm after 15 generations.

It interesting to monitor the solution process with a History plot of the objective function (see figure 4) where the generations can be easily highlighted. The moving average can be also plotted, tracking in some sense the convergence of the genetic algorithm to the solution.

The best solution found with NSGA-II after 15 generations is reported in table 1. This solution can be considered satisfactory from an engineering point of view; the difference between the experimental and simulated values of vertical displacements is actually sufficiently small for our purposes. The maximum is registered at the last load step and it is only 0.031 [mm].

Fig. 5: Deformed shape corresponding to the best solution found.

In figure 5 the final deformed shape corresponding to the best solution is plotted. It can be seen that the auto contact in the V-shaped profile of the rubber has been activated.

The principal strain component in the rubber ring ranges, more or less, between -0.7 (compression) and 0.7 (tension); the maximum shear strain has a maximum around 1.4. Therefore, the Mooney-Rivlin three terms model can be considered, in this context, as a good choice for an accurate description of the rubber behavior.

Conclusions
In this article a possible strategy to solve inverse problems with modeFRONTIER has been presented. As already mentioned, a similar procedure can be adopted for the identification of other interesting material properties, such as the viscoelastic ones.
However, it is worth mentioning that this procedure can be extended to the solution of other inverse problems arising in different engineering situations.
If the run time required by a single simulation makes the solution of the inverse problem prohibitive, the meta-modeling tools in modeFRONTIER could be used to drastically reduce the solution effort.