A Review of the Challenges and Limitations of Full‐Field Measurements Applied to Large Heterogeneous Deformations of Rubbers

Abstract:  This paper presents an overview of the use of full‐field measurement techniques, more precisely digital image correlation (DIC) and coupled DIC and infrared thermography, for the material and structure characterisation of rubber reported in the literature. Even though such techniques have increasingly been applied for approximately 30 years for moderate deformations in metal and composite materials, they are still under‐employed in the measurement of full kinematic and thermal fields in the case of large deformations undergone by rubber materials. To date, the applications addressed are crack propagation at both macroscopic and microscopic scales, model validation and constitutive parameter identification.


Introduction
Experimental mechanics applied to the field of rubber dates from the 19th century. The most cited works in the literature are those of Gough [1] and Joule [2], which dealt with the increase in temperature of rubber when stretched. At the end of the nineteenth century, mechanical tests were carried out to measure the change in volume of stretched rubber. The first record of this phenomenon dates back as far as the end of the 19th century in the works of Joule [3]. The author observed that the specific gravity of natural rubber (NR) decreased when it was stretched. But rubber mechanics really dates from the beginning of the 20th century, with a number of discoveries made at that time, such as the fact that a decrease in rubber stiffness is observed during the first mechanical cycles [4]. Later, this phenomenon was studied more precisely by Holt [5] and Mullins [6] and was finally referred to as 'the Mullins effect' (the reader can refer to Ref. [7] for further information on this effect). Then, while studying the nature of the stress-strain relationship for rubber containing different pigments in varying quantities, Schippel considered that possibly, when the rubber was sufficiently stretched, it might pull away from the particles of pigment along their axes of stress and cause vacua to be formed on either side of each particle, and a considerable increase in rubber body volume might therefore be observed [8]. Since the 1940s, the ability of rubber to resist fatigue loading has been investigated. Cadwell [9] and Beatty [10] were the pioneers in this field. Their work characterised the behaviour under fatigue loading of rubber samples and highlighted the beneficial effect of stress-induced crystallisation on life expectancy. The period between the 1950s and the 1980s was largely devoted to the development of theoretical works, and major results were obtained, including for instance the theory of rubber elasticity [11], the consequence of the stress-induced crystallisation on the behaviour of rubber [12] and the definition of tearing energy as the large strain counterpart of the elastic energy release rate [13].
The 1980s saw the appearance of numerical methods such as the finite element method [14] that allows mechanical engineers to predict local quantities, for example in the crack tip zone. Even if predictions using such methods can easily be compared to experimental data in the case of small deformations, to improve modelling, this was not possible in the case of large heterogeneous deformations. Indeed, due to the high deformability of rubber, the measurement of full mechanical fields such as the displacement field was not possible using classical experimental methods, typically gauges or a series of pyrometers. For this reason, the experimental measurement of local mechanical quantities was not really envisageable until the 1980s. By then, the improvements in image processing using numerical tools had enabled noncontact measurement techniques to become more and more popular in the experimental mechanics community. Currently, even though such optical measurement methods are commonly used to measure full fields at the surface of materials, they remain under-employed for the experimental mechanics of rubber. This is mainly explained by the fact that the large deformations undergone by such materials involve the development of specific methodologies for full-field measurements. This topic has recently been addressed in the literature, and it is now possible to characterise the mechanical response of rubber using digital image correlation (DIC) techniques. In the case of rubber material, these techniques have four main advantages: • the measurement is taken without any contact with the material, and therefore, it does not disturb the local displacement measurements; • the full kinematic field is measured. Thus, heterogeneities such as cracks are characterised in terms of heterogeneities in the kinematic field. The technique is therefore well adapted to the identification of mechanical properties from tests inducing heterogeneities in the deformation field; • a direct link can be established between the physical phenomena studied (their effect on the kinematic field) and continuum mechanics; • coupling full kinematic and thermal measurements leads to quantitative calorimetry and therefore to a tool for thermomechanical model validation and for studying the nature of various dissipative mechanisms.
A number of techniques are used to measure the full kinematic field. The nature of the measurements can be displacement or strain. Techniques used to measure displacement are numerous [15]: for instance, speckle [16]; speckle interferometry [17]; geometric moiré [18]; moiré interferometry [19]; holographic interferometry [20]; image correlation [21]; or the grid method [22]. Strain can be obtained by the numerical differentiation of the above displacement fields with suitable algorithms [23] or directly, for instance with shearography [24,25], speckle shearing photography [26] or by moiré fringe shifting [27]. Further information is given in Ref. [28].
In the case of rubber materials, techniques used have to conjugate many advantages: the measure of very large strains, the measure at very different scales, some facilities for the surface preparation and the possibility to use simultaneously an IR camera. Even if DIC technique is not well suitable to produce measurement for highly heterogeneous deformation due to a non-punctual measurement, it is the only technique that offers such a trade-off. This explains why the DIC technique seems to be the most suitable for characterising the large and heterogeneous deformation of rubber. DIC has proved over the last 20 years to be a very valuable tool for full-field displacement measurements. Even though DIC is classically used for moderate deformations of materials, typically inferior to 10%, this study will try to demonstrate that the DIC technique is also well adapted to rubber characterisation and is still under-employed for characterising and identifying rubber properties. The pioneering works in this field were those of Septanika et al. concerning the formalism of DIC in large deformations [29], Laraba-Abbes et al., applying DIC to the case of large homogeneous deformations [30], and Chevalier et al., who applied DIC to the case of large heterogeneous deformations [31]. Concerning the measurement of heterogeneous temperature fields, this is classically carried out using infrared cameras [32][33][34]. The features of this technique are not recalled here.
This study reviews studies of the mechanical and thermomechanical characterisations of rubber using the DIC technique, infrared thermography (IRT) and DIC coupled with IRT in the case of large heterogeneous deformation fields. First, mention is made of the special attention that must be paid on technical aspects for DIC and IR thermography. Second, applications dealing with crack growth in rubbers are reported at both microscopic and microscopic scales. Third, applications dealing with quantitative calorimetry applied to large and heterogeneous deformations of rubber are presented. Finally, full-field measurements during heterogeneous tests offer novel perspectives in terms of model validation and characterisation of dissipative mechanisms that are discussed in the last section of the study. Challenges and limitations are given for each application reported.

DIC and IR Thermography Applied to Rubber
Technical aspects for application of DIC to rubber: surface preparation and error evaluation The DIC technique requires sufficient black and white contrast to make the correlation between grey levels that are classically obtained by a speckle on the sample surface. Generally, it does not matter how the speckles are formed, which means that incoherent light can be used. For small deformations, the speckles are classically obtained by spraying a painting, white if the initial surface colour is black or inversely, the both colours successively if the initial surface does not provide enough contrast with only one painting. These techniques are not really suitable for large deformations undergone by rubbers. Typically, these techniques lead to a lack of correlation at large strains due to delamination between sprayed particles and rubber. To overcome this problem, several solutions are used to obtain a speckle with particles which are well distributed at the microscopic scale and consequently to obtain an optimised histogram in terms of greys levels. Several microscopic particles can be tested [35]: for instance, silica, clay, chalk, talcum powder, painting, but the main challenge is the choice of the method to deposit the particles: spraying, sprinkling, brushing, etc. Each of these solutions has to be checked with respect to the fact that DIC measurement errors strongly depend on the random speckle on the sample surface and on the fact that the good distribution of grey levels remains totally unchanged between two deformation states. It should be noted that an alternative method is reported in the literature for transparent silicone rubbers. It consists in generating speckle inside the rubber, on a plane between two rubber layers using the method developed in Ref. [36]. If no particles can be sprayed or deposited on the sample surface, typically when observing the interior view of a crack that opens and closes, DIC can nevertheless be performed using the natural speckle, and the lighting is the main adjustment parameter to improve grey level distribution and to reduce correlation errors. Consequently, significant preparation time is required. DIC error assessment is usually made in the case of homogeneous mechanical transformations, namely rigid body translations, planar rotations, or out-of-plane rigid body motions, which result in apparent essentially affine transformations of the 2D image [37][38][39][40][41][42]. Some authors use synthetic images that mimic real patterns to compute displacements in the Fourier [41] or space [38,39,43] domains or record a speckle image of an actual experiment and apply artificially imposed displacements [40,42]. For instance, in Ref. [42], the mean of the standard displacement uncertainties is founded to decrease as the subset size increases with a power law variation, showing that the displacement uncertainty and the spatial resolution are always the result of a compromise. For further information, the reader can refer to Ref. [44].
In the case of very large displacements of material points, the reference picture is generally updated and the measurement uncertainty will increase from analysis to analysis. This effect may induce too low signal-to-noise ratios. Conversely, if the pictures are not updated, the correlation will eventually fail and no results at all can be used. Therefore, the higher the picture increment, the lower the measurement uncertainty, provided the correlation results are trustworthy. This trade-off should be evaluated by assessing the performance of the DIC code. Due to the significant (artificial or natural) speckle variation, the value of the standard uncertainty level can significantly vary at each reference picture updating. For this purpose, the methods described previously could be applied. An alternative method consists in calculating the a priori performance of the correlation [40]. The central part of the Region of Interest is artificially moved by a displacement of 0.5 pixel in both directions using the shift/modulation property associated with Fourier transforms. This leads to the maximum standard uncertainty level for a given interrogation window size. As an updating procedure is used to evaluate large strain levels [45], this value is representative of each incremental measurement. This method was used in Ref. [46]. As an example, the standard displacement uncertainty was found to be equal to 0.035 pixel for an interrogation window size equal to 16 pixels.

Technical aspects on the surface preparation for IR thermography
Temperature measurement based on IR thermography is classically taken with a sample surface whose emissivity is close to 1. This level of emissivity is ensured by spraying a mat black painting on the sample surface. This preparation allows the emissivity to be at least superior to 0.93 [47][48][49][50]. An alternative technique can be founded in Ref. [51]. It consists in depositing carbon nanoparticles by means of candle combustion. However, for large deformations, either these techniques lead to a significant change in the surface emissivity at large deformations due to delamination between sprayed particles and rubber, or the coating is not chemically compatible with the rubber matrix, or the deposition technique is not suitable. For this reason, the full thermal field measurement is usually performed with carbon black-filled rubbers whose surface emissivity is close to 1. Nevertheless, the fact that the emissivity does not change with the strain level has to be checked for each type of filled rubber tested.
Surface preparation becomes more complicated when coupled full kinematic and thermal field measurement is envisaged on the same surface. Indeed, a satisfactory grey level distribution is needed for DIC while ensuring a surface emissivity close to 1. Thus, new coatings have been developed for this purpose: talcum powder combined with titanium, zinc, iron and chrome oxides and dissolved in ethanol [52]; white painting spread on a surface beforehand painted in mat black colour [53]. This type of surface preparation, which is mainly carried out on metallic materials, is not suitable for rubbers, and new measurement techniques have to be considered. They are precisely detailed in Calculation of heat sources in the case of large heterogeneous deformations.

DIC at the macroscopic scale
The hallmark of full-field measurements is to detect and characterise heterogeneity in the mechanical fields. For this reason, it is very suitable for providing information of importance on crack occurrence and growth in materials. Even though this kind of application is commonly addressed in the case of small deformations, the literature is very meagre concerning full-field measurements of large heterogeneous deformations undergone by rubbers. Such studies provide measurements of the kinematic field at the crack tip. Classically, results are discussed according to the evolution of the major strain with the distance to the crack tip [54][55][56]. Some of them account for the effect of the initial crack length [57] and that of the crack tip radius [58,59]. To the knowledge of the author, only that reported in Ref. [60] investigate the effect of temperature on the kinematic field at the crack tip to link it to the evolution of the microstructure. The authors use full-field measurements to characterise the effect of changes in the microstructure on the kinematic field at the crack tip. The high strain level encountered at the crack tip of crystallisable rubbers ensures that the material is partially crystallised. As the increase in temperature tends to melt the crystallites, it is interesting to characterise the effect of crystallite melting on the kinematic field. For this purpose, the authors use a rubber sheet, notched at its centre using a razor blade. When it is stretched, the crack opens and its shape becomes circular. A heat source is then applied to the crack lip via a heated metallic cylinder. Figure 1 gives information on the experiment and the results obtained with respect to the deformed sample geometry. This figure highlights that the local stretch ratio at the crack tip exceeds 1.64, i.e. the stretch ratio at which crystallisation begins. Contrarily, the stretch ratio obtained on either side of the crack tip is close to 1 and highlights the fact that these zones are relaxed. Moreover, the authors highlight the contribution of both internal energy and entropy variations to the changes in the kinematic field.
The optimisation of both lighting uniformity and speckle allows the kinematic field in the final deformed state to be obtained directly by correlating the grey levels between the undeformed and the final deformed states. This last remark is a key point of the measurement, because the reduction in the number of iterative steps avoids (in this study) or at least limits error accumulation -as discussed in Ref. [46] resulting from several correlations between the undeformed and the final deformed states. This highlights that the software used, and more particularly the algorithm used, plays a preponderant role on the quality of the final result. In this work, the resolution reached is 0.03 pixel, corresponding to 1.02 lm, and the spatial resolution attained (defined as the smallest distance between two independent points) is 10 pixels, corresponding to 341 lm. The software used for the correlation process is SeptD [61].
Given that large deformations are generated at the crack tip and that the scene size is large, the camera has to be located far from the sample surface to measure the displacement of the points located around the crack. Consequently, this leads to an increase in the resolution, but a decrease in the spatial resolution. A trade-off has therefore to be found to characterise accurately the high deformation gradient in the crack tip zone. It should be noted, finally, that for high deformation levels at the crack tip, most rubbers exhibit out-of-plane displacements that should be taken into account in the correlation process, for instance using profilometry [35], or using 3D DIC [62,63].
The DIC technique at the macroscopic scale offers various perspectives in terms of the validation of propagation laws and the prediction of kinematic fields at the crack tips. For instance, it provides interesting information on the distribution of the loading cases encountered in the crack tip vicinity. Their mapping during crack propagation could improve the understanding of crack growth mechanisms. Another perspective in this field is the use of the extended finite element method, firstly developed in the framework of finite element simulation [64] and recently applied to DIC measurement to account for discontinuities [65]. The method is based on a decomposition of the displacement field onto a regular finite element basis supported by a uniform mesh, enriched with suitable functions to describe accurately discontinuities.

IRT at the macroscopic scale
Similar to the DIC technique, IRT brings data of importance on heterogeneity in the kinematic field at the crack tip. Thus, studies mainly consist of mapping the full temperature field at the crack tip and discussing on the temperature gradient [56]. Mapping is obtained at a given loading, and the evolution of the temperature is never given for a particular point during the test. This is explained by the large displacement of the points on the observed surface. However, it allows us to identify regions within the sample that are undergoing accelerated damage under fatigue [66]. In this case, IRT is a suitable technique that enriches the analysis of crack propagation direction in a non-uniform stress field and as a function of both the angle of the initial cuts and the initiation of crack branching and catastrophic failure. Figure 2 gives an illustration of this kind of measurement. It should be noted that the thermal images are generally obtained with a pixel calibration procedure provided by the camera manufacturers and called the non-uniformity correction. Such calibration procedure considers the detector response to be linear with the radiation flux, and the measure of defective detectors (pixels) is replaced by that of detectors in their vicinity, which could be a problem in high thermal gradient zones. For these reasons, alternative procedures have been recently developed to account for a particular calibration law at each detector and to exclude the defective pixel from the thermal image [52,67].

DIC at the microscopic scale
Digital image correlation at smaller scales is difficult, because of the large displacement of the points at the  [60]. From the top left corner to the top right corner, the rubber specimen was previously cut with a razorblade before being stretched to put a heated metallic cylinder inside the crack. The full displacement field around the crack is obtained using the digital image correlation technique. The stretch ratio in the loading direction, denoted k in the diagram at the bottom of the figure, is then plotted versus the time at points A, B, C and D surface observed. Even if DIC at the microscopic and nanoscopic scales was recently performed at the surface and in internal planes of stretched rubbers [68,69], the displacements measured in these studies were low. At the nanoscale, for instance, the displacement applied to stretch the specimen was inferior to 100 nm and the scene size was 50 · 50 lm 2 . At the microscopic scale, the size of the scene was 500 · 400 lm 2 . Moreover, a speckle was applied to the surface observed. At the both scales, the displacement of the points of the considered plane was very low and significantly smaller than that at crack tips. The study of crack growth mechanisms requires the capacity to measure, at the microscopic scale, large displacements of the points on the observed surface. Moreover, the fact that during such experiments, the crack can be cyclically closed excludes the use of a speckle sprayed on the surface. This is why, to our knowledge, only one study addresses such a topic in the literature [70]. In this study, the authors used a long-distance microscope to observe the interior of the crack. The size of the scene was 2.7 · 2 mm 2 , and the authors measured the displacement of each zone during the stretching of the specimen. Thus, the observa-tion of crack propagation was enriched by the measurement of the relative displacement of the crack interior zones, as shown in Figure 3. The diagram in this figure gives the variation in the relative strain at different points of the crack tip during stretching. The term relative is used to define a quantity that is calculated from an intermediary strain state and not from the undeformed state. It is indeed impossible to use the undeformed state as the reference state because no measurement can be taken at the closed crack tip.
The observations made showed that from a highly stretched level, the crystallised crack tip does not deform any more (see relative strain at point B in the diagram) and resists crack growth. Coupled with further experiments, the authors identified the mechanism of fatigue crack growth in carbon black-filled NR submitted to severe loading conditions. Perspectives for such measurements are various; one is addressed in this study and deals with the fracture of carbon black agglomerates located in the crack path. Figure 4A shows a carbon black agglomerate which has failed at the fracture surface. The fracture surface of the agglomerate is parallel to that of the sample. Both are  [70] perpendicular to the cyclic loading direction. The question here is to define precisely how they fail, i.e. does the crack propagate through them or do they fail just before, due to stress concentration in the crack tip zone? To address this question, full kinematic field measurements can be taken on the crack tip interior using a CCD camera equipped with a long-distance microscope. Here, the natural speckle of the surface is used. After optimising the lighting, the resolution and the spatial resolution reach 0.03 pixel, corresponding to 0.06 lm and 16 pixels, corresponding to 31 lm, respectively. It should be noted that the spatial resolution attained corresponds to at least one-tenth of the average size of agglomerates (for example between 300 and 500 lm for N326 carbon blacks). This means that the displacement of the points in the vicinity of the agglomerate is precisely measured. The software used for the correlation process was Correli LMT [45]. During symmetrical stretching, the crack opens and propagates by generating surface from a tear-like line. Figure 4B presents the observed scene, corresponding to the stretched crack tip once the crack has propagated through the filler agglomerate. Figure 4D gives the full kinematic field (displacement in the loading direction in pixels) just before the agglomerate is revealed. The current image used for the correlation is given in Figure 4C. The reference image is that corresponding to the stretched crack tip just before detecting any strong heterogeneity in the agglomerate zone. As shown by this measurement, strong heterogeneity in the displacement field in the stretching direction is observed on both sides of a tear-like line. This is explained by the failure of the agglomerate due to the increase in stress concentration in this zone. This is therefore the evidence that in this case, the agglomerate fails before the crack propagates through it.
As a summary, full-field measurement techniques, and more particularly DIC, provide data of importance to improve the understanding and the modelling of crack growth in rubbers. At the macroscopic scale, they enable models to be validated by comparing measured and predicted heterogeneous fields. However, the fact that the measurement uncertainty increases from analysis to analysis (at each reference image updating) could be a limitation of DIC for large deformations of rubber, especially for the determination of the strain concentration at the crack tip. At the microscopic scale, full-field measurements enable the estimation of the relative strain of the zones that compose the crack tip. Thus, they enable physical mechanisms involved during crack growth to be identified. At this scale, the fact that the Zone of Interest is not perfectly plane and that out-of-plane displacements are generated during deformation leads to a significant misestimation of the kinematic field. This is the reason why DIC at the microscopic scale in rubber remains rather in the range of physical and phenomenological description and is not used quantitatively for identification steps.

Coupling Kinematic and Thermal Measurements: Towards Quantitative Calorimetry for Large Heterogeneous Deformations in Rubbers
To date, IRT has proved to be a relevant technique for studying engineering materials such as steel, aluminium and composites. Various studies previously carried out by Chrysochoos et al. have shown that heat sources produced by the material itself are more relevant than temperatures when analysing various phenomena such as Luder's bands [71], fatigue [72] or strain localisation [39]. A review of this research area is given in Ref. [73]. The main reason is that the temperature field is influenced by conduction as well as by heat exchanges with ambient air and grips, unlike local heat sources. The calculation of the latter requires tools dedicated to processing temperature maps provided by infrared cameras. Such tools have been developed in the case of small deformations, i.e. when current and initial configurations are not distinguishable, or in the case of large homogeneous deformations, and are no longer suitable for large heterogeneous deformations for which the infrared camera measures the temperature of points that move past each detector. It is therefore necessary to track each point during its large displacement to obtain the evolution of the temperature at each corresponding point. It seems that this issue has only seldom been addressed in the literature for rubbery materials. In Ref. [74] for instance, the authors propose two motion compensation techniques to account for this effect, but it seems that the strain amplitude is not significant, as only thermoelastic effects are studied, while calculations are carried out within the framework of large deformations. Nevertheless, studies carried out on steels and NiTi alloys have addressed the spatial synchronisation of kinematic and thermal fields [53,75]. The methods developed in such studies could easily be applied to large heterogeneous deformations of rubber.
In this section, the theoretical framework is briefly recalled before presenting the methodology used for calculating heat sources in the case of large heterogeneous deformations.

Theoretical background
The processing of the experimental 2D temperature map requires the use of successive hypotheses to write a bidimensional interpretation of the heat diffusion equation. First of all, in a thermomechanical framework [76], the local state axiom is assumed [77]. Any thermodynamic system out of equilibrium is considered as the sum of several homogeneous subsystems at equilibrium. The thermodynamic process is considered as a quasi-static phenomenon. The state of any material volume element is defined by N state variables: temperature T, a strain tensor denoted E and internal variables V 1 ,V 2 ,…,V N)2 such as plastic strain or volume fractions of some phases. The specific free energy potential is denoted W(T,E,V k ). Considering the first and second principles of thermodynamics and assuming that Fourier's law is used to model heat conduction, the heat diffusion equation can be written as follows: where q is the density, C E,V k is the specific heat at constant E and V k , K is the thermal conductivity tensor, and r is the external heat source. The right-hand side of Equation (1) represents the heat sources s produced by the material itself. It can be divided into two terms that differ in nature: • mechanical dissipation d 1 (or intrinsic dissipation): this positive quantity corresponds to the heat production due to some mechanical irreversibilities such as internal friction or crack occurrence and growth; • thermomechanical couplings: these correspond to the couplings between the temperature and the other state variables. The coupling between temperature and strain qTð@ 2 W=ð@T@EÞÞ _ E is the thermoelastic coupling [72,78]. If the temperature variation is governed by that of the internal energy, such as metallic materials, the thermoelastic coupling causes a temperature decrease (negative heat source) for a positive strain rate _ E and inversely. For rubbery materials, a competition takes place between the effects of the internal energy variation and the entropy variation with respect to the strain level, and a thermoelastic inversion can be observed [2]. At low strain levels for a positive strain rate _ E, the variation in temperature is governed by that of the internal energy and a temperature decrease (negative heat source) is observed. At higher strain levels for a positive strain rate, the variation in temperature is governed by that of the entropy and a temperature increase (positive heat source) is observed. Finally, it should be noted that coupling effects involving V k may create significant heat sources. For instance, if V k represents the volume fraction of a given phase k, the qTð@ 2 W=ð@T@V k ÞÞ _ V k quantity is a latent heat production [79].

Usual assumptions to calculate heat sources in the case of small deformations
The approach classically used to assess heat sources from the temperature fields obtained by infrared camera [80,81] is shortly described in this section.
Bidimensional formulation of the heat diffusion equation. It is first assumed that the problem is bidimensional. For the sake of simplicity, heat conduction will be considered as isotropic and Equation (1) can be rewritten as follows: where D is the Laplacian operator and k the conductivity coefficient which replaces the thermal conductivity tensor K in the case of isotropy. It is worth noting that this assumption is very strong, especially in the case of crystallisable rubbers. By integrating the heat diffusion Equation (2) over the sample thickness [82] and defining the mean thermal disequilibrium over the thickness between the sample and its surroundings by h ¼ T À T ref , where T ref is the temperature of the reference state, the following bidimensional formulation of the heat diffusion equation is obtained: where D 2D is the two-dimensional Laplace operator in the (x,y) plane. s can be considered as a time constant characterising the heat losses by convection and radiation between the sample surface and the surroundings. It is defined as follow: where e is the sample thickness, h is a convection coefficient. Note that s should be measured for each testing configuration (loading condition, machine used, environ-ment, etc). Moreover, it can also depend on the stretch ratio applied and the temperature range in which it is identified. The external heat source r is taken into account by monitoring T ref of reference state of the sample of the same geometry placed near the specimen or of the reference state of the sample itself.
To estimate the heat sources s produced by the material, the procedure consists in calculating the left-hand side of Equation (3) by processing the bidimensional temperature fields. It should be noted that the large deformations undergone by rubbers lead to discuss on the limitations of the thermophysical properties of the rubber and of the bidimensional formulation of the heat diffusion equation in terms of assumptions used. Indeed, emissivity is assumed to be constant and homogeneous at the sample surface during the deformation of rubber. This assumption has to be checked for each test performed. Moreover, strain induces anisotropy in the heat diffusion, due to chain orientation and chain crystallisation when possible. Also, the variation of the time constant with the strain level and with the range of temperature considered should be characterised. These are the next challenges to quantitatively study the thermomechanical couplings in such materials. Finally, due to the bad thermal conductivity of rubber, the sample thickness should be sufficiently small to assume the heat sources produced by the material itself to be constant through it.

Calculation of heat sources in the case of large heterogeneous deformations
The large deformations undergone by rubber require the dissociation of the initial and current configurations. If the deformation field is heterogeneous, the displacement of each point has to be tracked. For this purpose, both kinematic and thermal fields have to be measured, and several techniques can be employed. Two cameras (kinematic and thermal) can be placed one on either side of the sample, and a mirror symmetry can be performed between the scenes observed by the two cameras. If this is not possible, for instance due to the specific device used (see the example Ref. [60]), measurements could be taken on the same surface. In this case, the cameras are positioned perpendicularly and a dichroic filter-mirror allows each camera to observe the same scene [83]. Special attention must be paid on the speckle, which can disturb surface emissivity and therefore the temperature measurement. An alternative is the use of reflective spots with a different emissivity from that of the sample surface to build the displacement field from only the full thermal field using the same method as finite element approach. From a metrological point of view, the thermal resolution of the IR camera and the spatial resolution obtained have to be suitable to the measurement of the thermal gradient in the considered zone. Then, the main metrological difficulty is to measure the full kinematic field with a spatial resolution that is close to that of the IR camera. If such conditions are met, the temperature at a given point can be used quantitatively.
Concerning the calculation of heat sources, two strategies can be employed to differentiate spatially and temporally the experimental temperature field. Either the differentiation is performed in the initial configuration, and the lagrangian heat source field is calculated using the deformation gradient tensor to pull back the current temperature field, or it is performed in the current configuration to obtain the eulerian heat source field. In the former case, special attention must be paid to the temporal differentiation of the thermal field. Indeed, the temperature gradient and the displacement rate at each considered point come into play (contrary to the case of a homogeneous eulerian deformation field). Figure 5 shows results obtained using reflective spots at the surface of a notched elastomeric specimen submitted to sinusoidal mechanical cycles [84]. Seven zones are chosen to present the results in the vicinity of the crack. Figure 5 shows the location of the seven zones and their temperature variations. Processing the temperature fields by compensating for the movement of the points shows that temperature variations in zones 3, 4 and 5 are quite different from their counterparts obtained in zones 1, 2, 6 and 7. Zones 3, 4 and 5 are located near the crack tip, i.e. in a highly stretched area, and zones 1, 2, 6 and 7 are located in low-stretch zones. The temperature variations in zones 1 and 7, as well as zones 2 and 6, are quite similar. This can be explained by the quasi-symmetry in the displacement field on either side of the crack tip and by the fact that the zones are close to each other. Moreover, the temperature variation Figure 5: Temperature variations versus time for seven zones located around the crack tip of a rubber specimen submitted to cyclic mechanical loading. The thermal response is in phase with the displacement applied to the sample by the moving grip. The sinusoidal global stretch ratio varies between 1.03 and 1.67 [84] in zones 1 and 7 (i.e. the less stretched zones) clearly highlights the thermoelastic effects discussed previously (see Theoretical background). In fact, temperature does not change sinusoidally with time: during the loading, the temperature first decreases before increasing, and during unloading, temperature decreases before increasing at small strain levels, typically inferior to 7%. This phenomenon does not appear in zones 3, 4 and 5 in which the minimum strain level is higher. Indeed, the macroscopic stretch ratio variation (between 1.03 and 1.67) induces a microscopic stretch ratio variation (between 1.14 and 5) at the crack tip whose minimum value is superior to that at which thermoelastic inversion occurs. Figure 6 gives the heat source mapping calculated in the initial configuration.
Perspectives offered by quantitative calorimetry applied to large heterogeneous deformations of rubbers are promising and might improve significantly the knowledge of crack growth mechanisms. Moreover, it provides experimental data necessary to characterise the dissipative or coupled nature of the deformation mechanisms and to validate thermomechanical model predictions. Nevertheless, such analysis is complicated. First, from a metrological point of view, with respect to the temperature variation measured, the Narcissus effect significantly disturbs the measurement and must be subtracted from each image recorded by the infrared camera [84]. Second, the possible change in the emissivity of the surface observed has to be characterised for accurate heat source calculation. Similarly, constant s must be identified in the range of the temperature variation at the sample surface and its possible variation with stretching must also be characterised. Third, the assumption of the isotropy of the thermal conductivity is strong and may significantly disturb the calculation of heat sources, especially at high strain levels. Moreover, it should be noted that the thermal conductivity is generally not constant with respect to the temperature in vulcanised elastomers, and it should be necessary to account for its linear decrease with the temperature variation [85]. Finally, the resolution difference between the kinematic and thermal fields must be considered when trying to quantify precisely the heat source gradient in the crack tip zone.

Model Validation and Identification
Due to the large deformations undergone by rubbers, measurements that require contact with the material surface are excluded. This means that no experimental vali-dation of the heterogeneity induced in the kinematic and thermal fields was really possible until non-contact measurement methods were used. Full-field measurements consequently contribute greatly to advances in this field. The aim of this section is first to review the interest of fullfield measurements for model prediction validation and second to report the studies dealing with the identification of constitutive parameters from heterogeneous tests. A number of studies report the full kinematic field measurement in the case of homogeneous deformations (see for instance Ref. [30]). These studies, which tend to validate full-field measurement techniques as a relevant means of non-contact measurement, are not reported here. The present review focuses specifically on large heterogeneous deformations. Before discussing model validation and identification, a brief description of the specific behaviour of rubber is given.

Some specificities of the mechanical response of rubber
Rubber is a material with 'extraordinary physical properties', as Joule wrote in 1884 [3]. Most of these properties remain unexplained, this is why today rubber is still the object of vivid scientific interest. The main specificities of the mechanical behaviour of rubber can easily be highlighted by homogeneous uniaxial tensile test. Figure 7 shows the evolution of the nominal stress versus stretch ratio in carbon black-filled NR during five cycles performed at increasing strain levels. As shown in this figure, the behaviour of rubber is non-linear in terms of material behaviour. The material stiffness significantly decreases between the first and second cycles. The higher the strain applied, the higher the loss of stiffness. After a few cycles, the material can be considered as accommodated. No accommodation is observed for a strain level superior to that previously applied to the material. Moreover, during the second cycle, the stiffness increases for strain level close to the maximum strain level applied at the first cycle. Such a behaviour is further complicated by great sensitivity to the effects of loading rate, loading direction, temperature, environment, nature and amount of fillers, and method of manufacture. The occurrence of these numerous phenomena is clearly a difficulty for modelling, and currently, no  model accounts for all of them. Attempts to characterise the mechanical behaviour of rubber fall primarily into two categories: attempts based on statistical thermodynamics and those which take the phenomenological approach of treating the material as a continuum. The statistical thermodynamics approach [86] is based on observations that the rubber elastic force arises almost entirely from the decrease in entropy with increase in applied extension. This approach is considered to be only applicable to around 50% strain, as suggested in Ref. [87]. The majority of research works appear then to have focused on developments of the phenomenological approach, based on the observation of rubber under various conditions of homogeneous strain. Rubber is thus characterised by a strain energy density (SED) W, the stored strain energy per unit volume. This quantity is expressed as a polynomial function of strain invariants or directly in terms of principal stretch ratios (a review is given in Ref. [88]). Nevertheless, the complexity of the rubber behaviour makes it difficult to model the various phenomena involved in the deformation of rubber. Phenomena classically modelled are Mullins effect [89][90][91][92], anisotropy [93,94], viscosity [95][96][97][98], volume change [99], temperature effect [97,100] and stressinduced crystallisation [101,102]. The complexity in the formulation of most of these models is clearly a limitation for parameter identification and the inverse method to use. This is discussed in Identification from heterogeneous tests.

Validation of model predictions
To validate model predictions from heterogeneous tests, high gradients in the mechanical field are required, particularly due to the non-linear strain-stress relationship. Experimentally, several possibilities are offered to obtain such conditions. For instance, Meunier et al. have recently proposed to stretch a five-hole rubber sheet [103]. Two of the holes were linked together by cutting the rubber between them. Finite element simulations of the test were carried out using the four strain energy densities to be compared [104][105][106][107]. The predicted kinematic field and deformed geometry were compared with the experimental measurements and the results highlighted the limitations of certain constitutive equations were used to predict the deformation of highly stretched zones. A similar approach was adopted in Ref. [108] who tend to validate the prediction of their numerical simulation in the case of the bulge test [109]. The deformation field at the surface of the rubber sheet under such loading conditions is measured using an optical system [110], which measures tridimensional coordinates of circular markers previously painted on the specimen surface. The constitutive parameters used for the finite element simulation are fitted beforehand from homogeneous tests. Alshuth et al. have validated their model, which accounts for stress softening using full kinematic field measurements at the crack tip, and they compare the decrease in the major strain with the distance from the crack for both the first and fifth cycles to model predictions [57]. Finally, only one study deals with the validation of the prediction of heat sources at the crack tip using full temperature field measurements [66]. In this study, IR thermography was used to measure surface temperature profiles in filled NR and filled styrene butadiene rubber specimens during tensile fatigue cut growth tests. A linearised constitutive approach and finite element analysis were used to evaluate heat generation and associated transient temperature fields. Results obtained were compared with the IR measurements.

Identification from heterogeneous tests
The validation of model prediction involves fitting the constitutive parameters beforehand. Classically, identification is carried out from different homogenous tests, namely Uniaxial Tension (UT), Pure or Simple Shear (PS or SS) and Equibiaxial Tension (ET) tests. The latter test is often conducted using the heterogeneous bulge test, for which only the centre of the sample submitted to equibiaxial loading is considered [103,111,112]. The result of identification procedure is a trade-off between the parameter sets obtained separately with each homogeneous test. This identification procedure, which is widely employed, stems from the fact that in most cases, the strain energy densities chosen do not account for all the tridimensional phenomena involved in the deformation of rubber -viscosity, stress-induced crystallisation, manufacture-and deformation-induced anisotropy, stress-softening, thermoelastic inversion, to give some examples. In any case, the motivation for proposing identification from one heterogeneous test is to put several homogeneous tests together. Such heterogeneous tests result from the combination of sample geometry and loading conditions.
Four kinds of heterogeneous test are reported in the literature: (i) an inflated rubber sheet, namely the bulge test, which activates loading cases from ET to PS, theoretically, [109]; (ii) a cruciform rubber sheet, stretched (symmetrically or not) in the perpendicular directions. In this test, the branches of the sample are submitted to UT, the centre to ET and the intermediate regions to PS [46]; similarly, Drapier et al. have proposed a biaxial test based on the PS test [113]; (iii) a three-branched sample, stretched in the three branch directions, which leads to the same kind of induced heterogeneity as (ii). It is the only test performed with a conventional uniaxial tensile machine [114]; (iv) the tear test [115].
The analysis of test-induced heterogeneity is of the first order to obtain an efficient identification. Heterogeneity is driven by the combination of two criteria: (i) the number of different loading cases induced by the deformation of the specimen; (ii) for a given loading case, the nature of the loading level distribution. A discussion of this topic is proposed in Ref. [114]. Moreover, the effect of induced heterogeneity must also be discussed with respect to the considered identification procedure.
To identify the constitutive parameters from heterogeneous strain fields, an inverse problem has to be resolved. Parameter identification is performed within the framework of constitutive equations, which are a priori assumed to be relevant. The experimental data required are the displacement or strain fields, the sample geometry and the resulting forces applied to the sample. Three strategies are reported in the literature to address the inverse problem: updating finite element models, the virtual field method (VFM), and identification using DIC software.
The finite element method provides mechanical fields for almost any loading conditions, specimen geometry and constitutive equations, which explains its great popularity. It is moreover a tool for solving the inverse problem iteratively. A first calculation is carried out with an initial set of constitutive parameters. Nodal displacements are collected and compared to their experimental counterparts. The difference is quantified using an objective function to minimise iteratively with respect to constitutive parameters. This objective function is often the sum of the squared differences between the numerical and experimental data. The procedure stops when the objective function is lower than a given threshold value. Such procedures have been successfully simulated for the identification of constitutive parameters from a heterogenous tear test [113,115]. In the both studies, the strain energy density is a Rivlin series, i.e. a polynomial function of the strain invariants, with four parameters in Ref. [113] and three parameters in Ref. [115]. It should be noted that in Ref. [115], the viscosity is modelled by the Prony series and the Mullins effect is modelled with a damage function proposed in Ref. [89]. Nevertheless, the viscous and damage parameters are identified separately using homogeneous tests, only the hyperelastic parameters are identified from the heterogeneous test.
The virtual field method, proposed in Ref. [116], clearly departs from the finite element approach, as it is direct and does not require any iterative calculation in many cases. The basic idea is quite simple, as it consists of writing the global equilibrium of the tested specimen with the principle of virtual work. Introducing the constitutive equations and particular virtual fields leads to an equation where only the constitutive parameters are unknown, provided that the strain fields are measured, either directly or indirectly, by the numerical differentiation of measured displacement fields, for instance. The idea is then to write the principle of virtual work with different and independent virtual fields. This leads to a system of equations in which the constitutive parameters are unknown. If at least as many virtual fields as unknowns are chosen and if the constitutive parameters are written as polynomials of the strain, the system becomes linear and the unknowns are directly determined after the inversion of the system. This approach has recently been developed in the framework of hyperelasticity [46,111,114]. As an example in Ref. [46], parameter identification is carried out from a two directional tensile test performed on a cross-shaped elastomeric specimen. A virtual field has been constructed piecewise, using a set of twelve rectangular subregions shown in Figure 8A. In each of these subregions, the virtual displacement field is constructed using shape functions similar to those used in four-noded quadrangular finite elements [14]. As a result, the virtual displacement within each of these subregions depends on the nodal displacements. In the framework of the Mooney hyperelasticity [104], the optimisation of the set of two virtual fields provides two linear equations where the constitutive parameters C 1 and C 2 are unknown. Then, C 1 and C 2 are determined for different values of the global stretch ratio k g corresponding to each image stored, as shown in Figure 8B,C. Then, the final value of C 1 and C 2 is defined by their average values C 1 and C 2 as follows: where stdx ij (stdy ij , respectively) is the standard deviation of the difference between raw and smoothed u x (u y , respectively) distributions, i ¼ 1,…,2,j ¼ 1,…,n t , where n t is the number of images. The parameter set {C 1 ;C 2 } obtained is close to that obtained from the three homogeneous tests.
The main advantages of this method are first its directness in many cases and second the fact that the influence of the above problem posed by the boundary conditions can be avoided by choosing virtual fields in which only the resulting forces F virtually work [28,[117][118][119][120][121]. It should be noted that even if the virtual field method has only recently been applied to the identification of the elastic properties of rubber,it represents the largest proportion of the studies in this field, and numerous works are currently in progress in several laboratories. For instance, an interesting work has recently been published on the use of the virtual (A) (B) (C) Figure 8: Identification of constitutive parameters using the virtual field method (A) specimen (solid line) and zones defining the piecewise virtual field (dotted line) [46] (B) C 1 (C) C 2 field method for identifying the anisotropic and hyperelastic behaviour of in vitro human arteries from full-field measurements [122]. Identification can also be carried out using a DIC programme that evaluates the whole strain field at each step of the heterogeneous deformation, and the stress field is calculated from the estimated constitutive parameters. The forces in each direction can be processed and compared with the measured force in the loading direction. Constitutive parameters are then modified to fit the experimental data. This approach was used in Ref. [31] to identify the constitutive parameters of the hyperelastic model by Lambert-Diani and Rey [123], which is formulated using the strain invariants.
These identification procedures based on the resolution of an inverse problem, applied to large and heterogeneous deformations of rubber, offer a promising alternative to classical identification from several assumed homogeneous tests. However, much progress still needs to be made. For instance, constitutive equations used do not depend on time and temperature, and consequently, they do not account for the complexity of the physical phenomena involved in the deformation of rubber. Such considerations have to be addressed for competing with the classical identification procedure based on homogeneous tests. Nevertheless, the use of more complicated models could be a strong limitation for the resolution of the inverse problem using the virtual field method or the use of the DIC programme compared to the finite element method. Currently, only the exploratory study in Ref. [124] using the finite element method has envisaged the identification of the thermo-visco-hyperelastic behaviour of rubber from full kinematic and thermal field measurements. Finally, the use of quantitative data obtained with DIC for identification mainly concerns the macroscopic scale. As discussed in Technical aspects for application of DIC to rubber: surface preparation and error evaluation, the DIC error and the signal to noise ratios have to be considered for the identification of the constitutive parameters, especially if the reference picture is updated. DIC at the microscopic scale, or at the macroscopic scale but in high gradient zones, typically in the crack tip vicinity, is still in the range of physical and phenomenological description.

Conclusion
Full-field measurement techniques are still underemployed in the mechanics of rubber, although they offer new means of investigation in terms of mechanical characterisation and identification. The techniques reported in the literature are mainly DIC and IRT. At the macroscopic scale, the measurement of kinematic and thermal fields provides experimental data necessary to validate the prediction of thermomechanical models. The coupling between the two techniques offers interesting perspectives in terms of quantitative calorimetry. In particular, dissipative mechanisms such as viscosity, crack occurrence and growth, microstructure changes due to temperature variations can be characterised quantitatively. At the microscopic scale, the full kinematic field measurement provides information of importance to identify deformation and damage mechanisms involved in crack growth. It represents therefore complementary analysis to that carried out with scanning electron microscopy. Finally, the knowledge of full mechanical fields at the macroscopic scale allows us to identify parameters of constitutive equations by means of three main methods reported in the literature, namely updating finite element models, the virtual field method and identification from the DIC software itself. Currently, the constitutive equations used to model the behaviour of rubber for inverse identification from heterogeneous tests are not really realistic. The use of more sophisticated models is clearly a challenge for identification from heterogeneous tests to compete with the classical identification from several homogeneous tests.