Abstract
For studying the interaction of displacements, stresses, and acting forces for elastic and viscoelastic materials, it is of utmost importance to have a decent mathematical model available. Usually such a model consists of a coupled set of nonlinear differential equations together with appropriate boundary conditions. However, since the different material classes vary significantly with respect to their physical and mechanical behavior, the parameters which appear in these equations are unknown and therefore have to be determined before the equations can be used for further investigations or simulations. It is this very step which is addressed in this article where we consider elastic as well as viscoelastic material behavior. The idea is to compute the parameters as solutions of a minimization problem for Tikhonov functionals. Tikhonov regularization is a well-established solution technique for tackling inverse problems. On the one hand, it assures a computation that is stable with respect to noisy input data, and on the other hand, it involves desired a priori information on the solution. In this article we develop problem adapted Tikhonov functionals and prove that a Tikhonov regularization improves the accuracy especially when the underlying system is ill-conditioned.
1. Introduction
A material modeling process cannot be seen as a simple and single procedure, since it consists of multiple tasks. Usually, the beginning is the determination of the material behavior out of experimental data. In this step of a material modeling process, the data basis for the further work is established. It has to be decided if elastic, plastic, or viscose effects or combinations of these have to be included. The underlying structure of the model is established on this choice. The mathematical description has to follow physical principles, so that in a second step, the theoretical background from continuum mechanics has to be considered, typically leading to a model which consists of a coupled set of nonlinear differential equations. Subsequently, the outcoming material model has to be realized numerically. Based on these results, it is now possible to quantitatively compare the simulation to the experiments, which yields a parameter identification. In general, an inverse method is proposed, which determines the model parameters from the minimization of the error between model and experiment.
This last point, the parameter identification, is mainly examined in this contribution. Most of the attention is usually dedicated to the three firstly mentioned steps: experiment, theory, and numerics, whereby the final point of the identification process is often disregarded. Especially in some recent works, it comes up that the importance of the parameter finding increases more and more since the models try to include more material characteristics often resulting in a high number of model parameters. Moreover, some previous contributions have shown that a mechanical characterization only based on uniaxial data is not sufficient for a realistic description of three-dimensional, inhomogeneous problems. For details concerning the necessity of a multiaxial approach the reader is referred to the contributions of, for example, Baaser et al. [1–3], Johlitz and Diebels [4], and Seibert et al. [5]. In a model based only on uniaxial data the simulations for a parameter identification can be executed on ordinary geometries, such as simple cubes, whereby in a multiaxial description the complete specimen geometry and the resulting inhomogeneities have to be considered. This results in inverse calculations, where the detailed experimental conditions are reproduced in the numerics and subsequently compared to the simulation.
A similar effort is necessary when a geometrical structure has to be represented. This is the case when a complete geometry of an assembly is investigated [6, 7] or if a microstructure is examined, for example, for composites [8–10] or foams [11].
Both the increasing model complexity and the necessity of multiaxial approaches result in a more complex simulation which needs high effort to be solved. Finally, the computational costs are very high and lead to long-lasting computations. Regarding the parameter identification, many simulations have to be executed to find a matching parameter set. Hence, it is definitely recommended to treat the identification process much more attentively. An efficient parameter identification can reduce the duration for a finished modeling process with matching parameters a lot.
In a material modeling process it is usual to invest a high effort in the data acquisition in the experiments and the model description, whereby it is a common method to use only very simple optimisation algorithms for a parameter finding. Stochastic methods such as evolution strategies [12, 13] or the pattern search-algorithm [14, 15] as well as the method fminsearch based on the Nelder-Mead simplex algorithm [16] are often applied. The advantage of all of these methods is that no gradient information of the model is needed so that the effort in realizing and starting the parameter identification process is very low, whereby the performance of the algorithms is normally not satisfying. Concerning the realization it is only necessary to run the simulation and to define a fitness function which is able to compare the experimental and numerical results adequately. Due to the ill-posedness of the inverse problem, small errors in the input data of the identification may lead to large errors in the model parameters. This problem can be overcome by the introduction of appropriate regularization strategies. Regularization methods are stable solution strategies for inverse problems. They assure that the computed solution fits the given mathematical model and at the same time attenuate the noise that is contained in the measured data. Theory and application of regularization methods in Hilbert spaces are well-founded; we refer to the standard textbooks [17–19]. In the last decade the theory was extended to Banach spaces; see [20]. Amongst the most popular regularization tools are Tikhonov functionals [21]. These consist of two parts, a data fitting term and a penalty term. The latter is responsible for the stability and at the same time for the incorporation of any a priori information. That is why Tikhonov functionals might be an interesting tool to tackle parameter identification problems connected to constitutive equations of elastic and viscoelastic materials. This is the research hypothesis to be validated in the present article. We construct Tikhonov functionals that are adapted to our models and demonstrate the performance of Tikhonov regularization in different settings and condition numbers of the underlying system.
The applied methods and their advantages are shown with respect to a material model according to the work of Scheffer [22]. In order to focus on the process itself the data bases for the identification are not the original experimental values, where signal noise and statistics play an important role. Moreover, it is not finally sure whether the resulting parameters in the work of Scheffer [22] represent the optimal set since the model is of high complexity in order to describe a nonlinear viscoelastic material behavior depending on the deformation velocity. Hence, for this article the model with the already identified parameters is taken as the reference and the parameters are reidentified. As the main advantage the exact values of the parameters are known and it can be proven if the proposed methods find this exact solution of the problem.
Organization of the Article. In Section 2 the theoretical aspects concerning the material description and the resulting constitutive equations are discussed. In Section 3 the problem of parameter identification as an inverse problem is addressed. Section 4 verifies the usage of the proposed methods by showing theoretical results for simplified material models as well as by reconstructing given material parameters. The paper is concluded with a summary and an outline of future research in Section 5.
2. Mathematical Model
A new approach for identifying material parameters is discussed here. There are many viscoelastic materials, and different models exist to characterize them fully. To test the approach, the following model is used: It is inspired by real life experiments and is based on a rheological model, where a spring is positioned in parallel to Maxwell elements. Each of the Maxwell elements consists of a spring connected in series to a dashpot; compare, Figure 1 (Figure 1 is reprinted by permission from Springer Customer Service Centre GmbH: [23]).
The single spring represents the equilibrium stiffness, whereas the Maxwell elements describe the viscoelastic effects. Even the 3-parameter model consisting of a spring and one Maxwell element in parallel shows the main viscoelastic effects, namely, relaxation and creep. The mathematical structure of this rheological model forms the basis for the nonlinear and three-dimensional formulation used in the following investigations.
The basic elasticity is given as the equilibrium response and serves as the basis for the further material model. A basic elasticity curve for elastomers often shows stiffening for larger deformation values resulting in the typical S-shaped stress-strain curve. This shape of the stress-strain correlation can be monitored by a Yeoh model [24] with the free energy function:with the first invariant in the third order. Here is the left Cauchy-Green deformation tensor and is the deformation gradient. The higher-order approach can monitor the S-shaped upturn of the stress-strain correlation for bigger strains. Figure 2 shows a typical stress-strain diagram for a carbon black filled rubber sample for the basic elasticity as investigated in [23].
To cover the different time ranges in the relaxation behavior, four Maxwell elements are used. The viscoelastic behavior can be captured by the easiest neo-Hookean approach for the springs in the Maxwell elements. This leads to the free energy contribution:The total free energy is given by the sum of (1) and (2):In the equations, describes the deformation in the single equilibrium spring and represents the elastic deformation in the th Maxwell element. The underlying multiplicative split of the deformation gradient was well established in finite strain theories of plasticity or viscoelasticity; compare, for example, [25]. The Cauchy stress is obtained from the free energy functions according to the Clausius-Planck inequality for isothermal processes with respect to the argumentation of Coleman and Noll [26] leading towith as the Lagrangian multiplicator representing the incompressibility of the material. Evaluating (4) for the chosen free energy function results inThe time-dependent behavior of the Maxwell elements is described by the evolution equationsfor the inelastic deformation in the dashpot of the th Maxwell element; see [27]. are the material time rates of the inelastic right Cauchy-Green deformation tensors . The parameters in the evolution equation stand for the so-called relaxation time. They have to be chosen according to observation from relaxation tests. Usually different relaxation processes take place in a black carbon filled rubber. While the fastest process is triggered by the loading time , the slowest process is governed by the total time of observation . For simplicity the relaxation times of the model are equally distributed with this time interval on a logarithmic scale. Figure 3 shows the theoretical strain that a material experiences during the relaxation experiments and its stress answer.
(a)
(b)
In the experiments, the stretch values , which are the ratio of current length to initial length , are steadily increased until the maximum values are reached and are then held there. The time when this maximum is reached is denoted by . While in an idealised relaxation test the loading is applied instantaneously, in real experiments is mainly governed by the maximum speed of the experimental device. Furthermore, the stop time of the experiments is denoted by . In the stress answer of the material an increase is noticed until where the maximum is reached. From then on, the stress tends towards the basic stress value while the deformation does not change. This process is called relaxation.
The relaxation times were chosen in such a way that the physical restrictions in the experiments are preserved. On one side the maximal device velocity that is used in the experiments gives a lower bound for the smallest observable relaxation time. On the other side the time span of the experiment determines the longest relaxation time. In order to cover a broad interval, the relaxation times are set in a range of multiple powers of ten as seen in Table 1. It also shows the identified material parameters from the aforementioned experiments performed on a black carbon filled rubber [23]. They are mechanically meaningful and can be used to test the here proposed reconstruction methods. Due to the complexity of the model the identification of its parameters is far from being trivial. In general, an inverse problem is formulated in such a way that the parameters are chosen by minimization of the error between data and experiments.
3. Regularization Methods
If we neglect modeling errors, then the exact parameters and those which have been predicted by solving the inverse identification problem, coincide. This leads to the equationwhere are the measurements at different time steps , the approximative model, and the unknown (exact) material parameters. Due to inevitable errors in the experiments and the fact that the model is unable to describe the real process exactly, the residue can only be approximated byWith the goal in mind to minimize the residues with respect to the model parameters, the and are unified into vectors and , respectively, and one obtains the minimization problemDue to the nature of the presented constitutive model (5), the equations are simplified to linear ones with respect to with a common form. This leads for values of the stretch toFor different values of the stretch , which explicitly depends on time, we obtain the stress as a linear combination of the coefficients determined from the model equations and the material parameters . The coefficients are given by (5) and are united into a vector :For rate-dependent material behavior, the stress and the coefficients in have explicit additional dependency on the inelastic stretch .
In the following we do not investigate data measured in real experiments but for simplicity we treat the so-called problem of parameter reidentification. Therefore, the used data are not derived from actual experiments but are obtained by the stress-strain relationship given by the constitutive equations and a set of known parameters; compare Table 1. To further simulate errors in the measurement process, a perturbation by white additive Gaussian noise is added to the simulated stress vector, which leads towhere is the simulated stress vector and the perturbed one. Here, measures the deviation between the two vectors and is called noise level. From here on out we use as (simulated) input data.
The vectors are collected to build the rows of a single matrix :Thus we have with . With the vector containing material parameters and containing stress values , one obtains the following minimization problem:with .
The minimization problem is differentiable due to the Euclidean norm and can therefore be solved by setting the first derivative equal to zero as a necessary condition. By defining a functionand its derivativesit can be shown that the -matrix is symmetric and positive definite, if has full rank, so that the extremum of is an absolute minimum. Setting leads to the so-called normal equationThis equation has a unique solution, if . This can be realized by using a sufficient number of samples. The solution of (17) is then given aswhere is defined as the pseudoinverseThe difficulties that may arise when computing the inverse of a matrix should be mentioned here. Explicitly, one should consider the condition of the problem. It tells how much the obtained solution deviates from the real one if the input data is perturbed. In general, a problem is called well-conditioned when small errors in the input data only generate small errors in the solution. On the other hand, an ill-conditioned problem gives rise to big errors in the solution for small perturbation of the data. For linear problems, the condition number of a matrix can be calculated with the help of the singular values bywith and being the largest and smallest singular value, respectively.
To compensate for the measurement errors which possibly lead to the bad condition of , we propose a regularization scheme that is based on the Tikhonov-Phillips method; see, for example, [18]. A penalty term is added to the minimization problem (14) such that the norm of the solution vector is to be as small as possible. This leads to the minimization problemwith the regularization parameter . By differentiating with respect to , we see that any solution has to satisfy the shifted normal equationSince the matrix is invertible, there exists a unique solution of (21) for all . The nature of the given problem encourages an additional penalty term. The second parameter of the basic parameters, that is, , has to be negative due to the characteristic shape of the stress-strain relation. Positive values of are therefore penalized such that (21) becomeswhere denotes a further regularization parameter. The optimality condition for the minimizing is then given aswith for and otherwise. The multiparameter setting as in (23) is closely related to elastic-net regularization which originates from statistics, see [28], and has been analyzed in [29].
4. Virtual Experiments
In the virtual experiments the three aforementioned regularization methods are used. Method 1 is the unregularized normal equation (17), method 2 is the regularized Tikhonov equation (22), and method 3 is given by (24), which adds an extra penalty term to guarantee the negativity of the basic elastic parameter . Furthermore we divide the experiments into three categories depending on the characteristics of the reconstructed parameters. At first, set 1 with the basic elastic parameters , , and is reconstructed. Subsequently the viscoelastic parameter set 2 with , , , and is computed using the elastic parameters as a priori information. Lastly, all parameters from the previous sets are put together in set 3 and reconstructed at the same time. As long as viscoelastic parameters are involved, we additionally consider a simplified model to work out additional dependence on the used relaxation times.
To compare the different reconstruction methods, an error measurement is introduced. The reconstruction error is the normalized difference between the vector with the identified parameters from Table 1 and the vector with the newly reconstructed parameters with this approach in the Euclidean norm; that is,In our tests we choose a noise level varying in the range of to display the results of our methods depending on the noise level. The value corresponds to exact (simulated) data. We perform 100 realizations for every value of which are averaged afterwards in order to reduce fluctuations in due to the added random noise. Note that the type of material being tested may also affect the noise level.
4.1. Reconstruction of Parameter Set 1: Reconstructing the Basic Elasticity Parameters
At first the basic elastic parameters are reconstructed. As our investigations are based on uniaxial tension tests, the tensorial constitutive equation is reduced to the scalar-valued stress-response as a function of the stretch:By evaluating (26) times, one obtains the linear system of equationswith , , and the perturbed stress vector .
In a first step, the condition of matrix has to be considered. Figure 4 shows the condition for different values of the maximum stretch . For a doubled probe length , the condition number decays rapidly if more data points are used until it is almost constant with . In the other cases, similar observations can be made and therefore is constant for further considerations. Since and the total time of observation are fixed, it is quite natural that the condition of saturates for growing at a certain level. The reason is that for a fixed number of columns and fixed observation time we do not expect that or as . Moreover, we have that , independently of , and thus there are at most singular values which are not identical .
We discuss now the regularization for two different values of the stretch, namely, and . The results are illustrated in Figure 5. A regularization for a material with double the length and does not seem to be necessary. Even the unregularized reconstruction gives reasonable low errors . For small derivations from the original stress vector , the results get even worse. For the smaller strain value the advantages of the regularization become more pronounced. The reconstruction error barely reaches a value of in method 2, whereas we face larger errors for the unregularized reconstruction. Additionally the advantages of the second regularization term with in method 3 can be shown. For the relative error may be small in method 2, but by a closer look it is revealed that 50 percent of the values of are positive and therefore the reconstruction is not satisfying. With the additional regularization term, the percentage of positive values can be reduced to 4 percent and the relative error is even further reduced. For the number of positive values for is also reduced by half in method 3.
(a)
(b)
The quantitative reconstruction of the parameters is shown in Table 2 for . For both choices of a sufficient reconstruction is possible with a worse result in the smaller interval up to .
4.2. Reconstruction of Parameter Set 2: Incorporating the Viscoelastic Parameters
The model is extended by incorporating viscoelastic parameters. The constitutive equation has the following form:As the identified basic elastic parameters are kept fixed in this step we only consider the viscoelastic parameter reconstruction and the data is solely produced by .
In the viscoelastic case the data points have to be distributed logarithmically on the interval , such that more weight is put onto the loading process of the deformation than the relaxation process.
4.2.1. Simplified Parameter Reconstruction
Due to the more complex material behavior, the importance of the relaxation times of the different Maxwell elements has to be considered first. To test the influence of the relaxation times on the reconstruction, simulated data with a simplified material model in mind have to be constructed. Only two Maxwell elements with relaxation times s and s are used. With this simplified model the influence of the different parameters can be studied.
The condition only differs slightly if the values for , or the number of data points is changed. On the other side, different values for have a big effect on the condition as can be seen in Table 3. Physically, for a given the parameter stays in direct relationship to the deformation velocity of the material. Slower deformation rates correspond to greater values of , whereas small values belong to faster deformations.
Table 4 shows the reconstruction of the simplified model if the parameters , are chosen as 1 MPa. A bigger value for the choice of the regularization parameter penalizes the norm of too much and should therefore be chosen accordingly. For further investigations a value of is chosen. It can be recognized that the reconstructed value of parameter is within a tolerance as soon as is smaller than .
Further observations can be made if the parameters are chosen similar to the relaxation times in decades. Table 5 shows the reconstruction for MPa and MPa and vice versa. In addition to the prior observations, it is also noted that the reconstruction of the parameters belonging to small relaxation times is not as good as the ones belonging to the bigger relaxation times.
4.2.2. Original Parameter Reconstruction
For the reconstruction of the exact parameters of the original model, the condition can be observed in Figure 6 for different values of . The number of samples is chosen as in these experiments. The actual reconstruction is shown in Figure 7 with different values for . According to earlier observation the regularization is not necessary for s. The error for the unregularized solution is as low as the error for the regularized one. The reconstruction error is small in either method though. A small means at the same time a relatively fast deformation rate. These rates cannot be obtained in current experiments. Therefore a more realistic value of s is also tested. In this case, regularization leads to better results again. The regularized reconstruction error never exceeds , whereas the unregularized one increases for bigger noise levels.
(a) s
(b) s
The quantitative reconstruction in Table 6 shows the characteristics of the simplified tests. For s, the parameter is reconstructed well, whereas the others show overshooting values. The smaller the relaxation time, the bigger the error gets. For s, which is closer to real experiments, the parameters get reconstructed appropriately, but a negative value for is observed. This behavior is according to the choice of as .
4.3. Reconstruction of Parameter Set 3: All-at-Once Computation of the Parameters
At last, the question arises if it is necessary to do the time-consuming identification of basic elasticity parameters beforehand. In case of highly filled rubbers the relaxation of the specimen can be observed for very long times. Therefore, the experiments are really time-consuming if the basic elasticity should be observed; see, for example [22]. In this case the sequential identification presented above is not possible. A simultaneous identification by means of our methods is proposed in the next step. For the following experiments the data vector consists of values of (28). The parameters and s are fixed. The regularization parameters are chosen as and .
4.3.1. Simplified Parameter Reconstruction
In the first step the simplified model from Section 4.2.1 is extended. Next to the parameters and for the viscoelasticity, the basic elasticity parameters , , and can be chosen freely. The condition of the -matrix is shown in Table 7. For values of as s or s, we still obtain a reasonable condition whereas we observe an exponentially increasing condition number as . To restrict the workload on the experimentator, the number of samples is chosen as . The values of the exact parameters have been chosen as in Table 8.
Table 9 shows reconstructions for different values of . It has to be distinguished between two types of parameters. On the one hand the viscoelastic parameters show the same patterns with respect to and the relaxation times as before. On the other hand an insufficient reconstruction of the basic parameters for a small choice of can be observed. Therefore the choice for builds a contradiction. One explanation can be found in the unusually high condition in the range of .
Furthermore, similar observations as before can be made if one of the viscoelastic parameters gets multiplied by 10, as seen in Table 10. By the Tikhonov regularization a simultaneous reconstruction of all parameters is performed since the minimization of (23) averages the norm of the penalty term. This leads to over- and underestimations in the reconstructed solution . A possible remedy could be a separate regularization with containing the basic elasticity parameters and containing the viscoelasticity parameters.
4.3.2. Original Parameter Reconstruction
In a last step the quality of the reconstruction of the actual parameters is tested. The range in which the condition lies does not change. For s the values lie in the range of going afterwards down to for larger values of . Due to this ill-posed problem a reasonable reconstruction without regularization is impossible for . With our two methods of regularization, however, sufficient results with an error measurement are reached.
Figure 8 shows the reconstruction of the parameter for s. Without regularization the error is increasing for higher values of . The regularized solutions never exceed . The importance of the regularization parameter is again shown. Even though the error is minimally larger for method 3, only 10 percent of the values for are positive whereas we have 53 percent positive values in method 2.
Table 11 shows the actual reconstructed values for two values of . We see the contradicting influence of . For a small value we reconstruct the viscoelastic parameters well but suffer a massive undershoot for the basic parameter . For the bigger value of the basic elasticity parameters can be reconstructed while the viscoelastic ones are slightly undershot.
Remark 1. Knowing the ground truth as in our simulations the regularization parameters and can easily be adjusted by trial and error. Of course, this is no longer possible in case of real, measured data. The most popular parameter choice rule in this case is Morozov’s discrepancy principle, which works as follows: Choosing null sequences and as well as a tolerance parameter , we define bywhere solves (24) with parameters , . Finally we set , . In this way it is guaranteed that the noise in the residuum is of the same order as the noise in the data. Heuristic methods for choosing , depending only on the measured data , which is of relevance if the noise level is not known and cannot be estimated, are the L-curve criterion or generalized cross validation. We refer the interested reader to the standard textbook [18] for more details on parameter choice rules.
5. Conclusion
In this article we introduced a method for identifying material parameters from constitutive equations by means of a Tikhonov regularization. An extra penalty term which is directly inspired from the stress-strain relationship for the chosen material is added to expect reasonable results. Overall, this method ensures a fast and convenient way to identify both the basic elastic and the viscoelastic parameters at the same time without a time-consuming separated identification. Maybe better reconstruction results can be achieved by using more sophisticated penalty terms. We furthermore mention some shortcomings that are intrinsic to Tikhonov regularization. The first one is the fact that in general yielding a much worse condition number for , the matrix that appears in systems (22) and (24). A second one must be seen in the determination of appropriate regularization parameters and in case of real data which might be time-consuming and data-dependent. Furthermore, better results can be achieved by correlating the relaxation times with the actual deformation process.
Finally we note that the methodology is also applicable to alternative constitutive equations like plasticity models. This will be subject of future research.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.