# Modélisation du comportement élasto-plastique des tôles par un modèle compartimenté

Laurent Tabourot, Pascale Balland, Laurent Bizet, Ahmed Maati, Virginie Pouzols

► To cite this version:

Laurent Tabourot, Pascale Balland, Laurent Bizet, Ahmed Maati, Virginie Pouzols. Modélisation du comportement élasto-plastique des tôles par un modèle compartimenté. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717040

HAL Id: hal-01717040

https://hal.science/hal-01717040v1

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Modélisation du comportement élasto-plastique des tôles par un modèle compartimenté

Laurent TABOUROT<fnref n="1" />*, Pascale BALLAND<fnref n="2" />, Laurent BIZET<fnref n="3" />, Ahmed MAATI<fnref n="4" />, Virginie Pouzols<fnref n="5" />

**Résumé** — Cet article a pour but de montrer que les incertitudes liées à la caractérisation du comportement élastique peuvent entraîner des variations significatives des résultats de simulations numériques. Il faut alors certainement introduire des phénomènes relativement fins pour disposer d'une description qui restreigne cette incertitude. Il y alors un intérêt à accroître le nombre de variables d'état des modèles pour raffiner la description du comportement élasto-plastique des tôles minces métalliques. Un modèle de ce type est appliqué, dans un premier temps, au cas d'une simulation de pliage d'un vé montrant que l'on dispose d'une marge de manœuvre pour mieux simuler des opérations de mise en forme de tôles présentant un grand retour élastique.

**Mots clés** — modélisation, comportement, procédé, hétérogénéités

## 1 Introduction

Les géométries d'un poinçon et d'une matrice utilisés pour la mise en forme des métaux en feuilles ne peuvent pas se définir directement à partir de la forme souhaitée du composant à produire. Ceci est dû à l'écart qui se produit systématiquement entre la pièce formée et l'outillage après le retrait de ce dernier en raison de l'inévitable retour élastique du métal déformé. Pour obtenir des pièces par formage conformes aux exigences du dessin de définition, il faut alors définir des formes d'outillage qui anticipent ce retour élastique. Comme il n'y a généralement pas de lien direct entre une modification sur l'outillage et sa conséquence sur le retour élastique, la mise au point par ajustements itératifs de la forme des outils est extrêmement laborieuse et ce d'autant plus qu'il est impossible de rajouter de la matière une fois celle-ci enlevée. La simulation numérique qui permet d'effectuer virtuellement l'opération de mise en forme et donc de multiplier à moindre frais les essais de géométrie apparaît comme une alternative intéressante pour la mise au point des outillages. La qualité de la modélisation du comportement mécanique du matériau est alors un élément déterminant pour la précision des simulations numériques.

La modélisation du comportement mécanique des matériaux métalliques nécessite de déterminer la relation entre le tenseur des contraintes appliquées et le tenseur des déformations produites. Elle suppose une partition tranchée de la déformation entre une déformation dite réversible (ou élastique) et une déformation permanente (ou plastique). La déformation élastique est imputée aux distorsions du réseau cristallin qui disparaissent lorsque la charge appliquée est supprimée. Les déformations permanentes résultent des réarrangements atomiques non réversibles après disparition du chargement. Ceux-ci sont induits par le mouvement de défauts linéaires (les dislocations) en réponse au chargement

<footnote n="1">1 Université de Savoie, Laboratoire Symme, France, laurent.tabourot@univ-savoie.fr</footnote>

<footnote n="2">2 Université de Savoie, Laboratoire Symme, France, pascale.balland@univ-savoie.fr</footnote>

<footnote n="3">3 Université de Savoie, Laboratoire Symme, France, laurent.bizet@univ-savoie.fr</footnote>

<footnote n="4">4 Université Amar Telidji, Laboratoire de génie des procédés, Algérie, taleb662003@yahoo.fr</footnote>

<footnote n="5">5 RAyce Eurl, France, Virginie.Pouzols@rayce-lab.com</footnote>

<footnote>* Auteur correspondant</footnote>
appliqué. Comme les déformations élastiques sont supposées dépendre linéairement des contraintes alors que les déformations plastiques apparaissent comme étant fortement non-linéaires et dépendantes du trajet de déformation (c'est-à-dire plutôt compliquées à appréhender), beaucoup de travaux de recherche se focalisent principalement sur le développement des modèles de comportement plastique. Avec la géométrie et les conditions aux limites, c'est le comportement plastique qui détermine principalement la distribution des déformations dans la pièce mise en forme. La modélisation du comportement élastique, avec son caractère linéaire a priori "simple" est donc, quant à elle, reléguée au second plan. On se contente généralement de supposer que le comportement élastique est linéaire, isotrope, indépendant des paramètres de mise en forme (température, déformation) et défini par les valeurs du module d'Young E et du coefficient de Poisson v la plupart du temps prises dans la littérature.

En dépit des efforts réguliers pour l'améliorer les modèles de plasticité et les codes de simulation numérique, la précision des résultats reste pourtant en deçà d'un seuil qui permettrait l'exploitation de ceux-ci à des fins pratiques d'optimisation des outillages dans un cadre industriel. Il devient alors raisonnable de rechercher d'autres pistes d'amélioration. L'étude de la sensibilité des résultats au comportement élastique du matériau est une voie qu'il convient d'explorer puisqu'il joue un rôle essentiel dans le retour élastique et donc pour sa prédiction.

L'objectif du travail rapporté dans cet article est donc de faire le point sur l'influence de la modélisation du comportement élastique sur le retour élastique du matériau. La caractérisation expérimentale du comportement élastique est effectuée sur un acier de type C68. On recense alors tout ce qui est source d'écart par rapport au comportement typique indiqué plus haut. On détermine alors l'incidence de ces écarts au moyen d'une simulation numérique sur la forme finale d'un vé obtenue par pliage.

On montre tout d'abord en recourant à trois méthodes de caractérisation qu'il est légitime de s'intéresser aux effets d'une possible dispersion des valeurs élastiques du matériau non déformé. Ces dernières sont clairement dépendantes de la méthode d'identification (et sans aucun doute de l'expérimentateur). Par un plan d'expériences, on détermine l'incidence de variations de ± 10% des paramètres élastiques sur la forme finale du cas de pliage considéré.

Le même travail est effectué en prenant en compte la décroissance d'un paramètre élastique au cours de la déformation observée expérimentalement.

Dans une dernière partie, on s'attache à analyser l'impact d'une remise en cause d'un comportement idéal linéaire purement élastique en-dessous d'une contrainte seuil. Il existe en effet plusieurs phénomènes susceptibles d'altérer la linéarité du comportement. Le recours à une modélisation permettant d'intégrer l'hétérogénéité du comportement est alors une première étape vers la prise en compte de phénomènes fins liés aux déformations réversibles.

# 2 Incertitude des constantes élastiques et conséquences

## 2.1 Détermination expérimentale du module d'Young

Il s'agit de montrer ici qu'il est raisonnable de considérer une dispersion des valeurs élastiques. Pour cela, le module d'Young est mesuré par trois méthodes différentes sur une tôle d'acier C68.

La première méthode consiste à déterminer le module d'Young comme la pente de la partie linéaire de la courbe contrainte déformation obtenue lors d'une traction sur une éprouvette. De façon à éviter les effets de désalignement de mors, de rattrapage de jeux sur la machine de traction, plusieurs charges - décharges (fig. 1) à faible déformation sont réalisées et la valeur de la pente moyenne est retenue comme valeur du module d'Young E (tableau 1).
<figure><img src="image_2.png" /><figcaption>Fig. 1: cycles de charge - décharge appliqués à une éprouvette en acier C68.</figcaption></figure>

<table><caption>Table 1 : valeur du module d'Young déterminée à partir de la réponse d'une éprouvette C68 à des cycles de charge - décharge en traction.</caption><thead><tr><td>Module d'Young (GPa)</td><td>Écart-type (GPa)</td></tr></thead><tbody><tr><td>212</td><td>1,56</td></tr></tbody></table>

La deuxième méthode consiste à mesurer la réponse vibratoire d'un échantillon du matériau suite à un choc et d'identifier les fréquences de résonance du matériau [1,2]. Cette technique permet d'évaluer non seulement le module d'Young mais aussi le coefficient de friction, le module de cisaillement et le coefficient de Poisson.

Pour réaliser une mesure, l'échantillon est donc soutenu par un système de fixation en 2 points et frappé par un marteau. Le signal sonore émis par l'échantillon en vibration est enregistré par un microphone et analysé pour en déduire les paramètres de comportement élastique. Cette étude expérimentale a été menée avec un dispositif de type RFDA Basic de la société IMCE. Une série de 10 essais a été réalisée sur un échantillon de C68 rectangulaire de dimensions 34,55×10,03×1,51 mm³.

<table><caption>Table 2: valeurs des paramètres élastiques mesurés à partir de l'analyse d'un signal obtenu par choc sur une éprouvette en acier C68.</caption><thead><tr><td>E (GPa)</td><td>Écart-type de E (GPa)</td><td>v</td><td>Ecart type de v</td></tr></thead><tbody><tr><td>217</td><td>0,02</td><td>0,28</td><td>0,0002</td></tr></tbody></table>

Le tableau 2 présente les paramètres obtenus après cette série d'essais. Une excellente répétabilité est constatée. La valeur du module d'Young est supérieure à celle déterminée par essai mécanique.

Enfin, la troisième méthode recourt à une mesure du temps de propagation d'une onde ultrasonore. Les vitesses de propagation des ondes dans un milieu élastique sont directement fonctions des paramètres du modèle de comportement élastique et de la masse volumique ρ du matériau. L'essai qui en découle consiste à mesurer ces vitesses de propagation afin d'accéder directement aux paramètres du comportement élastique. Une campagne de mesure a été menée avec un dispositif de mesure de temps de propagations des ultrasons de type Olympus 38DL Plus. À cet appareil se raccorde un transducteur distinct pour la mesure des vitesses de propagation des ondes longitudinales ou transverses. Les fréquences des ondes sont de 10 MHz pour les longitudinales et de 5 MHz pour les transverses. Des transducteurs à contact sont utilisés. Ils jouent à la fois le rôle d'émetteur et de récepteur. Le temps de vol est mesuré après réflexion de l'onde émettrice sur la face intérieure de l'échantillon.

Une série d'essais a été réalisée sur un échantillon d'acier C68. Le tableau 3 présente les résultats obtenus. Le faible écart-type indique une bonne répétabilité des mesures. La valeur du module est sensiblement inférieure à celle déterminée par essai mécanique. Cet écart peut
être dû aux incertitudes de mesure de la géométrie de l'échantillon mais aussi à sa faible épaisseur qui ne correspond pas aux conditions optimales d'utilisation pour ce type d'appareil (plusieurs millimètres). En effet, plus la distance parcourue par les ondes est importante, meilleure serala précision de la mesure. La caractérisation de tôles nous contraignant à travailler avec de faibles épaisseurs (de l'ordre de 1,5 mm), cette méthode n'est sans doute pas la mieux adaptée au cas du matériau caractérisé ici.

<table><caption>Table 3 : valeur du module d'Young obtenu par le temps de propagation d'une onde ultrasonore sur une éprouvette en acier C68.</caption><thead><tr><td>Module d'Young (GPa)</td><td>Écart-type (GPa)</td></tr></thead><tbody><tr><td>195</td><td>1,18</td></tr></tbody></table>

Les différences entre les valeurs obtenues par ces trois méthodes définissent l'incertitude qu'il convient d'appliquer aux paramètres élastiques. L'origine de cette différence doit pouvoir s'expliquer mais on souhaite surtout déterminer l'effet induit par cette variation expérimentale des valeurs. Cette étude n'a porté que sur la valeur du module d'Young mais on peut supposer que la détermination du coefficient de Poisson serait aussi problématique. Toutefois, cette étude expérimentale n'est pas traitée ici. La valeur de v est ici fixée par sa valeur couramment utilisée (0,3).

## 2.2 Incidence de la variation de la valeur initiale de E et v sur le retour élastique

Cette incidence est montrée par une simulation numérique d'un pliage avec un angle important pour obtenir des retours élastiques importants. Cette simulation est présentée sur la figure 2. Elle est effectuée avec le logiciel Abaqus.

L'éprouvette est en acier C68 (épaisseur e=0,5mm, longueur L=70mm, largeur l=10mm). Elle est maillée avec des éléments de type C3D8R. Le maillage comporte 19320 éléments dont 6 éléments dans l'épaisseur. Le poinçon possède un diamètre de 25 mm. La méthode d'intégration est implicite. Le poinçon est d'abord descendu sur une profondeur de 6 mm puis remonté. On a choisi volontairement d'imposer une flèche pas trop importante, de façon à avoir un grand rayon de pliage, ce qui conduit à des plus grands retours élastiques. On mesure alors l'angle d'ouverture θ après retrait du poinçon. Le poinçon, ainsi que la matrice en vé, sont modélisés par des solides analytiques rigides. Dans les simulations suivantes, la gestion du contact se fera avec un coefficient de frottement nul.

L'angle obtenu pour toutes les simulations, lorsque le poinçon est au maximum de sa course, est de 141,25°. La flèche sera définie par la différence de position entre l'état initial et l'état après retour élastique pour l'élément situé à L/2 et 1/2 de la surface supérieure de l'éprouvette.

<figure><img src="image_3.png" /><figcaption>Fig. 2: présentation de la simulation numérique d'un pliage en vé. Figure gauche au maximum de la course du poinçon. Figure droite après retrait du poinçon. L'angle d'ouverture θ du vé est mesuré dans cette dernière configuration.</figcaption></figure>
Pour ce travail, on suppose que le matériau est initialement isotrope et que l'écrouissage l'est également. La courbe de référence est celle obtenue par un essai de traction sur l'acier C68. L'incidence des valeurs initiales de E et v est testée en imposant des variations de ± 10% comme montré dans le tableau 4.

<table><caption>Table 4 : valeurs de (E, v) et (θ, flèche) correspondantes calculées par simulation.</caption><thead><tr><td>E (GPa)</td><td>v</td><td>θ (°)</td><td>Flèche (mm)</td></tr></thead><tbody><tr><td>195,3</td><td>0,27</td><td>167,3</td><td>-2,31</td></tr><tr><td>195,3</td><td>0,3</td><td>167,2</td><td>-2,33</td></tr><tr><td>195,3</td><td>0,33</td><td>167,1</td><td>-2,35</td></tr><tr><td>217</td><td>0,27</td><td>165,6</td><td>-2,63</td></tr><tr><td>217</td><td>0,3</td><td>165,5</td><td>-2,65</td></tr><tr><td>217</td><td>0,33</td><td>165,3</td><td>-2,67</td></tr><tr><td>238,7</td><td>0,27</td><td>164,1</td><td>-2,89</td></tr><tr><td>238,7</td><td>0,3</td><td>164,0</td><td>-2,92</td></tr><tr><td>238,7</td><td>0,33</td><td>163,9</td><td>-2,94</td></tr></tbody></table>

Dans ce même tableau, les valeurs de l'angle d'ouverture ainsi que de la flèche résiduelle sont reportées. Il convient de noter que proportionnellement une dispersion sur E a plus d'effet en pliage que la même dispersion relative sur v. La variation est de 0,15° par pourcent de variation sur E alors qu'elle n'est que de 0,01° par pourcent de variation sur v. La variation qui résulte des dispersions des trois mesures expérimentales ci-dessus serait donc de 2° ce qui est au-delà de la précision souhaitée pour ce genre d'opérations dans l'industrie. Il est également intéressant de constater que plus le module d'Young est important, plus le retour élastique est faible ce qui est bien conforme aux observations usuelles.

# 3 Dépendance du module élastique avec la déformation

## 3.1 Détermination expérimentale du module d'Young en fonction de la déformation

Des essais de traction avec cycles de charge-décharge ont été réalisés afin de déterminer le degré de dépendance du module d'Young avec la déformation plastique de l'acier C68. Les modules sont déterminés sur des éprouvettes ayant déjà subi une prédéformation de 2% à 20% avec un incrément de 2% [3]. La figure 3 présente l'évolution du module d'Young en fonction de la déformation plastique. Ces résultats montrent une diminution non négligeable du module d'Young de l'ordre de 10% pour une déformation plastique de 20% avec un module initial de 212 GPa.

<figure><img src="image_4.png" /><figcaption>Fig. 3: variation du module d'Young en fonction de la déformation plastique sur une éprouvette d'acier C68.</figcaption></figure>
## 3.2 Incidence de la variation de E au cours de la déformation

On recourt à l'utilisation d'une routine spécifique pour l'affection de E en fonction de la valeur de la déformation plastique. Les résultats sont donnés dans le tableau 5.

<table><caption>Table 5 : valeurs de l'angle d'ouverture et de la flèche résiduelle du vé obtenues par simulation numérique avec la variation de E au cours de la déformation plastique.</caption><thead><tr><td>θ (°)</td><td>Flèche (mm)</td></tr></thead><tbody><tr><td>166</td><td>-2,58</td></tr></tbody></table>

On constate que la valeur de la flèche correspond par extrapolation dans le tableau 4 à une valeur de E de l'ordre de 212 GPa qui est la valeur initiale. En fait la déformation imposée au niveau de l'éprouvette est relativement faible (déformation plastique équivalente maximum de 2,3%) et donc n'a pas d'effet sur le pliage tel qu'il est opéré dans le cadre de ce travail. Toutefois, il est vraisemblable qu'avec un petit rayon de pliage, des grandes déformations auraient été générées et une différence aurait été détectée entre les calculs avec ou sans prise en compte de la dépendance de E avec la déformation plastique. La valeur du module d'Young étant plus faible à grande déformation, on aurait tendance à avoir des retours élastiques plus importants dans la simulation avec la prise en compte de cette dépendance.

# 4 Effet d'une non linéarité du comportement élastique

## 4.1 Origine physique du comportement élastique

Il existe différentes raisons pour lesquelles le matériau n'a pas une raideur constante lorsqu'il est sollicité dans le domaine des déformations purement élastiques. La première de ces raison est que le matériau n'est pas indemne de défauts. Les dislocations sont largement distribuées à l'intérieur du matériau. Les segments de dislocations dont la longueur est variable se courbent sous l'effet de la contrainte appliquée. Ce phénomène est réversible pour les dislocations dont le rayon de courbure reste au-dessus d'une valeur limite. Cette pseudo-élasticité est d'une raideur moindre que l'élasticité du réseau. Il existe également un autre mécanisme corolaire du précédent. En effet, les dislocations courbées au-delà du rayon critique contribuent à la déformation plastique mais requièrent des contraintes moins importantes à déformation équivalente qu'une fraction de matériau déformée de façon parfaitement élastique. Par ailleurs, ces phénomènes sont à l'origine d'un stockage d'énergie élastique car des déformations élastiques supérieures à la déformation élastique moyenne du réseau sont requises pour accommoder localement l'écart entre la déformation de la zone qui s'est déformation plastiquement et la zone en déformation élastique pure. Ils sont également fondamentalement à l'origine de l'écrouissage cinématique. Ces phénomènes ne peuvent pas être pris en compte par une modélisation classique.

Une modélisation proposée récemment [4] permet de prendre en compte partiellement la microplasticité locale tout en permettant de considérer une valeur "idéale" du module d'Young. Cette modélisation considère fondamentalement le matériau comme étant hétérogène. Concrètement, chaque élément d'un maillage éléments finis se voit alors attribuer un modèle de comportement spécifique. Un cas assez simple consiste à définir alors des ensembles de lois de comportement bilinéaires dont l'écrouissage est constant et le deuxième niveau d'écrouissage identique. Seule la contrainte seuil est distribuée selon une distribution de type Rayleigh. Celle-ci a donc la forme mathématique suivante :

$$f\left(\frac{x}{\theta}\right) = \frac{2x}{\theta^2}e^{-\frac{x^2}{\theta^2}}$$
où $\theta$ est le paramètre à fixer qui pilote à la fois la moyenne et l'écart type de la distribution et $x$ est la variable variant de 0 à l'infini pour définir la distribution.

Sur la figure 4, la distribution de type Rayleigh ainsi que la courbe de référence décrivant un écrouissage bilinéaire distribuée en chaque point d'intégration du code de calcul par éléments finis sont présentées.

<figure><img src="image_5.png" /><figcaption>Fig. 4 : exemple schématique d'une distribution de type Rayleigh et courbes d'écrouissage correspondantes, toutes définies par une même pente et une même valeur de saturation.</figcaption></figure>

L'identification est menée sur un essai de traction. La moyenne des contraintes seuil initiales vaut 552,8 MPa, la valeur de la contrainte de saturations 934,2 MPa et la pente d'écrouissage n, 364,2 MPa. Dans ces conditions, il est intéressant de considérer le comportement du matériau dans la zone "élastique". On constate alors que le comportement du matériau n'est plus tout à fait linéaire : l'activation successive de différentes zones du matériau conduit à un affaiblissement progressif de la zone élastique jusqu'à l'entrée en plastification de l'ensemble des éléments, conduisant à l'apparition d'un écrouissage prononcé aux alentours de ce qui peut apparaître comme une contrainte seuil (400 MPa). La courbe contrainte déformation est toutefois correctement restituée au niveau de l'écrouissage.

<figure><img src="image_6.png" /><figcaption>Fig. 5 : allure de la courbe contrainte / déformation dans le domaine élastique.</figcaption></figure>

Cette distribution a été appliquée aux éléments composants le maillage de l'éprouvette du pliage en vé. Le même chargement a été appliqué. Les résultats sont donnés dans le tableau 6.

<table><caption>Table 6 : valeurs de l'angle d'ouverture et de la flèche résiduelle du vé obtenus par simulation numérique avec prise en compte des hétérogénéités du matériau.</caption><thead><tr><th>E (GPa)</th><th>ν</th><th>σ<sub>moy</sub> (MPa)</th><th>σ<sub>sat</sub> (MPa)</th><th>n</th><th>θ (°)</th><th>Flèche (mm)</th></tr></thead><tbody><tr><td>217</td><td>0,3</td><td>552,8</td><td>934,2</td><td>364,2</td><td>162,9</td><td>-2,92</td></tr></tbody></table>
On constate que la prise en compte des hétérogénéités donne une variation de 3° sur l'angle. Bien que le module d'Young apparent soit affaibli par rapport à sa valeur initiale de 217 GPa.

# 5. Conclusion

Dans cet article, l'effet du comportement élastique sur les simulations numériques est analysé. On montre que les incertitudes sur les paramètres élastiques suffisent à expliquer les différences notables entre la simulation numérique et les résultats expérimentaux : une variation de E qui entre dans les incertitudes de détermination entraine une variation d'angle supérieure aux tolérances attendues sur une pièce en pliage. Il faut affiner la description du comportement élastique pour disposer de marges permettant de régler son effet au plus juste sur le retour élastique.

Une façon originale d'intégrer l'effet de la microplasticité au moyen d'une distribution des contraintes seuils est proposée. La modélisation proposée est physiquement fondée. Son intérêt principal est qu'elle permet d'utiliser sans complexité supplémentaire les modèles usuels des codes EF et de restituer des comportements fins de comportement élastique, incluant l'effet de microplasticité et d'écrouissage cinématique. Ce raffinement autour du comportement élastique est important pour parvenir à prédire le comportement de pièces lors d'opérations impliquant de forts retours élastiques comme l'étirage-pliage par exemple.

# Références

[1] B.A. Latella, S.R. Humphries, Young's modulus of a 2.25Cr-1Mo steel at elevated temperature, Scripta Materialia. 51 (2004) 635-639.

[2] G. Roebben, B. Bollen, A. Brebels, J. Van Humbeeck, O. Van der Biest, Impulse excitation apparatus to measure resonant frequencies, elastic moduli, and internal friction at room and high temperature, Review of Scientific Instruments. 68 (1997) 4511.

[3] M. Vautrot, P. Balland, O.S. Hopperstad, L. Tabourot, J. Raujol-Veille, F. Toussaint, Characterization of high carbon steel C68 at elevated temperatures and different strain rates, Technische Mechanik. 32, 2-5, (2012) 577-586.

[4] L. Tabourot, P. Balland, J. Raujol-Veillé, M. Vautrot, C. Déprés, F. Toussaint, Compartmentalized model for the mechanical behavior of titanium, Key Engineering Materials. 504-506 (2012) 673-678.