# Approches local/global non-intrusives et optimisation robuste

Léa Karaouni, Ludovic Chamoin, Olivier Allix, S. Guinard

► To cite this version:

Léa Karaouni, Ludovic Chamoin, Olivier Allix, S. Guinard. Approches local/global non-intrusives et optimisation robuste. 26ème Congrès Français de Mécanique 2025, Aug 2025, Metz, France. hal-05108326

HAL Id: hal-05108326

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

Submitted on 11 Jun 2025

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 Attribution 4.0 International License
# Approches local/global non-intrusives et optimisation robuste

L. Karaouni<sup>a</sup>, L. Chamoin<sup>a</sup>, O. Allix<sup>a</sup>, S. Guinard<sup>b</sup>

a. Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS- Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France.

b. Airbus Group Innovations, Bât. Campus engineering, BP90112, Blagnac Cedex 31703, France.

lea.karaouni@ens-paris-saclay.fr, ludovic.chamoin@ens-paris-saclay.fr, olivier.allix@ens-paris-saclay.fr, stephane.guinard@airbus.com

## Résumé :

Les approches de submodeling sont très prisées par les industriels car elles permettent de réanalyser à faible coût, au sein de logiciels commerciaux, des zones locales jugées critiques en raffinant les maillages et la description de la géométrie ou même en adaptant les modèles de comportement. Une des limites majeures de cette méthode est d'être une approche purement descendante dont les résultats ne prennent pas en compte les changements globaux induits par les modifications locales qui peuvent être très diverses. Pour pallier cet inconvénient, les approches local/global non-intrusives ont été proposées afin de corriger les défauts du submodeling tout en restant compatible avec l'usage de codes industriels. Ce travail s'intéresse à l'extension de ces méthodes de couplage non-intrusives au cas de la conception robuste. Dans le cadre des approches stochastiques, on vise à propager les incertitudes locales au niveau global à l'aide de techniques non-intrusives afin de caractériser les distributions probabilistes des quantités globales induites par les incertitudes locales. Il s'agit en particulier de pouvoir proposer une optimisation globale robuste tenant compte des variabilités locales. Pour cela, nous utilisons l'algorithme génétique élitiste et non dominé évolué (NSGA-II) qui permet une grande exhaustivité dans la recherche de solutions optimales tout en préservant une complexité calculatoire acceptable. Cette approche permet notamment de gagner en efficacité lors de la phase de propagation d'incertitudes à travers le modèle multiéchelle. Dans ce travail, diverses méthodes stochastiques de résolution seront analysées et comparées (méthode par perturbation, chaos polynomial, Monte-Carlo...) sur des cas-tests académiques de difficulté variée.

## Abstract :

Submodeling approaches are widely used in industry, as they enable low-cost reanalysis within commercial software of local regions considered critical, by refining the mesh, enhancing the geometric description, or even adapting the constitutive models. A major limitation of this method is that it is purely top-down, and its results do not account for global changes induced by local modifications. To address this issue, non-intrusive local/global approaches have been proposed. This work focuses on extending such non-intrusive coupling methods to the context of robust design. In stochastic approaches, the goal is to propagate local uncertainties to the global level using non-intrusive techniques, in order to characterize the probabilistic distributions of global quantities that arise from local variabilities. In
particular, the aim is to enable robust global optimization that takes these local uncertainties into account. To this end, we use the elitist non-dominated sorting genetic algorithm (NSGA-II), which allows for an extensive exploration of the optimal solution space while maintaining reasonable computational cost. This approach significantly improves efficiency in the uncertainty propagation phase within the multi-scale model. In this work, various stochastic methods (perturbation method, polynomial chaos, Monte-Carlo, etc.) are analyzed and compared on academic test cases.

# Mots clefs : Conception robuste, Quantification et propagation d'incertitudes, Couplage local/global, Modélisation multiéchelle, Eléments finis stochastiques

## 1 Introduction

Les structures complexes industrielles sont généralement représentées par un système de modèles interconnectés plutôt que par un seul modèle unifié, avec la nécessité d'une communication entre les différentes échelles pour garantir la cohérence et l'exactitude du processus de simulation. Dans ce contexte, l'analyse des modèles en fonction de la maturité globale du projet est généralement menée de manière séquentielle et partitionnée (Figure 1). La conception globale est d'abord définie, puis l'architecture de second niveau est précisée, et ainsi de suite. Cependant, cette approche descendante peut facilement conduire à des conceptions conservatrices ou même non conservatives. En effet, l'impact potentiel des détails locaux (par exemple, modifications des types de connexions, positionnement des composants, etc.) sur l'architecture globale est généralement compensé par l'utilisation de coefficients de sécurité. Ce papier constitue une première tentative de proposer une approche complémentaire ascendante qui prend en compte, dès les phases préliminaires de conception, l'impact potentiel des changements locaux sur l'architecture globale. C'est la philosophie préconisée dans la notion de conception robuste, qui vise à gagner en robustesse vis-à-vis des variabilités et incertitudes [1–3].

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Processus classique de conception et flux associés</figcaption></figure>

Néanmoins, la description particulière du problème nécessite d'adapter la chaîne de simulation afin de permettre la communication entre les différentes échelles. Pour cela, l'approche usuelle est celle du *submodeling* car elle permet de réanalyser à faible coût, au sein de logiciels commerciaux, des zones
locales jugées critiques en raffinant les maillages et la description de la géométrie ou même en adaptant
 les modèles de comportement. Une des limitations majeures de cette méthode est d'être une approche
 purement descendante dont les résultats ne prennent pas en compte les changements globaux induits par
 les modifications locales et donc conduire à des erreurs importantes aux différentes échelles et à des
décisions de conception non conservatives. C'est pourquoi nous nous appuyons ici sur des approches
local/global non-intrusives proposées dans [4]. Ces approches permettent d'améliorer un modèle glo-
bal grossier sans altérer ses paramètres intrinsèques (matériaux, topologie, discrétisation, techniques
de calcul, nature du modèle, etc.) grâce à des patchs locaux contenant les complexités et variabilités.
Les données d'interface sont ensuite échangées de manière itérative entre les différents modèles. Cette
stratégie numérique constitue ainsi un outil qui corrige les inconvénients de la méthode de *submode-
ling* tout en restant compatible avec l'utilisation des codes industriels. Ce travail étudie l'extension des
méthodes de couplage non-intrusives à la conception robuste. Un premier pas dans ce sens concerne la
propagation d'incertitudes locales et globales dans le cadre du *submodeling* [5]. Notamment, un défi-
clé dans ce contexte est la propagation des incertitudes à travers les différentes échelles afin d'assurer
des prédictions de performance fiables. Pour cela, un cadre probabiliste est proposé pour propager les
incertitudes exploitant pleinement la nature non-intrusive de la méthodologie de couplage. L'expansion
en chaos polynomial et la méthode de perturbation sont envisagées comme des alternatives efficaces à
la simulation de Monte-Carlo. Dans ce travail, différentes techniques de propagation des incertitudes
seront comparées et testées sur un exemple illustratif en 1D particulièrement sensible aux perturbations
géométriques.

## 2 Couplage local/global non-intrusif

La stratégie multiéchelle proposée concerne les méthodes de couplage local/global non-intrusives [4]. Ces approches permettent d'enrichir un modèle global par des modifications locales sans altérer le modèle d'origine. Ces modifications peuvent concerner le comportement des matériaux, la topologie des composants étudiés, les techniques de calcul ou encore la nature même des modèles. Pour simplifier les notations, le principe du couplage est explicité dans le cas d'un couplage entre un tenseur de Hooke hérogène local et un autre homogénéisé global. Aussi, il est possible de considérer que les phénomènes d'intérêt sont localisés dans l'espace. L'idée est donc de décomposer le domaine initial Ω en deux parties sans recouvrement (Figure 2) :

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Décomposition de la géométrie en une zone locale et une zone globale</figcaption></figure>

- une zone locale Ω_l ⊂ Ω qui englobe le support du phénomène d'intérêt à analyser. Par souci de simplicité, le domaine local Ω_l est supposé strictement inclus dans Ω. Dans cette zone, un modèle basé sur la loi de comportement initiale complexe est conservé.

- la zone complémentaire Ω_c = Ω \ Ω_l dans laquelle un modèle plus grossier est pris en compte. Il est défini en remplaçant le comportement initial par un comportement linéaire élastique homogène.
L'interface entre les deux domaines $Ω_l$ et $Ω_c$ est notée $\Gamma$, et les modèles global et local sont supposés être géométriquement et cinématiquement compatibles sur cette interface $\Gamma$. Plutôt que de définir la solution comme la somme de la réponse globale et d'un terme d'enrichissement local, le problème local est intégré dans la formulation globale en utilisant une approche de substitution par morceaux. La solution locale remplace donc directement la partie correspondante de la réponse globale. Pour cela, l'approche de couplage local/global non-intrusive modifie le problème global en définissant le support de sa solution sur l'ensemble du domaine $Ω$, tout en prolongeant virtuellement le comportement homogène élastique linéaire sur $Ω_l$, sans modifier le problème local. Afin de respecter les exigences de non-intrusivité de la méthodologie, une méthode de Newton modifiée sur les quantités d'interface est généralement employée pour résoudre le problème. À partir d'une solution élastique globale, chaque itération se déroule comme suit :

1. **Analyse locale** : Une analyse détaillée est effectuée sur le modèle local, en imposant les déplacements globaux courants comme conditions aux limites.

2. **Calcul du résidu** : Le résidu est exprimé sous forme d'une force d'interface mesurant la violation des équations de couplage sur $\Gamma$. La convergence de l'algorithme itératif est vérifiée à cette étape. Si la norme du résidu est suffisamment faible, les itérations s'arrêtent.

3. **Correction globale** : Le résidu est réinjecté dans le modèle global sous forme d'une force interne supplémentaire sur l'interface, déformant ainsi le modèle global comme si les détails locaux étaient présents. La loi de comportement élastique linéaire et les conditions aux limites du modèle global sont préservées. Des techniques optionnelles d'accélération de convergence peuvent être appliquées à cette étape, puis la solution globale est mise à jour et le processus est répété jusqu'à convergence. Il est également intéressant de noter que la première itération de la méthode de couplage non-intrusive, sans réinjection de la force de rééquilibrage d'interface, correspond à l'approche classique de *submodeling*.

Cette formulation garantit la compatibilité cinématique entre les déplacements globaux et locaux ainsi que l'équilibre des efforts sur l'interface $\Gamma$. Cette classe d'algorithmes présente deux propriétés importantes. Premièrement, si l'algorithme converge alors sa limite est la solution de référence, i.e. la solution obtenue en effectuant réellement la substitution des modèles et en exécutant une analyse complète [4]. Ainsi, la technique de couplage non-intrusive est exacte et fiable : l'erreur peut être estimée à travers la norme des résidus et réduite autant que nécessaire. Deuxièmement, la méthode est en effet non-intrusive car les solveurs fonctionnent comme des boîtes noires et les données du modèle ne sont jamais modifiées car elle nécessite uniquement l'envoi et la lecture des quantités d'interface. Cela permet d'introduire facilement des lois de comportement complexes, des détails géométriques et des non-linéarités dans le modèle local, et donc d'effectuer une analyse structurale exacte sans avoir à modifier le modèle global. Ainsi, les techniques de couplage non-intrusives s'avèrent très utiles pour prendre en compte diverses modifications locales d'une pièce et sans nécessiter de remaillage. Notamment, dans notre cas, elles présentent un intérêt considérable pour la prise en compte de nombreuses variabilités locales.

## 3 Adaptation à la conception robuste

Les concepts et fondements de la conception robuste ont été développés dans les années 50 par Taguchi afin d'améliorer la productivité en ingénierie. L'objectif de la conception robuste est d'optimiser la performance moyenne du système et de minimiser sa variabilité (i.e. de diminuer la sensibilité) résultant des incertitudes dans les propriétés matérielles ou géométriques, sans en éliminer les causes. Or, les
problèmes rencontrés en ingénierie mécanique sont souvent des problèmes multiobjectif. Dans le cadre de la conception mécanique, on cherche à satisfaire des critères de performance antagonistes de sorte qu'il n'existe pas une solution unique permettant d'optimiser simultanément tous les critères de performance. L'idée est d'optimiser à la fois la valeur moyenne de la fonction de performance et sa variance en respectant les contraintes du problème considéré. Ainsi, le problème de conception robuste s'écrit sous la forme d'un problème d'optimisation multiobjectif et non déterministe.

$$ \text{Minimiser : } \mu_f(\mathbf{x}, \mathbf{s}) = \mathbb{E}_{\mathbf{z}} \left[ \max_{\mathbf{s} \in V} (\mathbf{f}(\mathbf{s}, \mathbf{z}, \mathbf{x})) \right], \sigma_f(\mathbf{x}, \mathbf{s}) = \mathbb{E}_{\mathbf{z}} \left[ \max_{\mathbf{s} \in V} \left( (\mathbf{f}(\mathbf{s}, \mathbf{z}, \mathbf{x}) - \mu_f(\mathbf{x}, \mathbf{s}))^2 \right) \right] \quad (1) $$

$$ \text{Sous les contraintes : } g_j(\mathbf{s}, \mathbf{z}, \mathbf{x}) \le 0, j \in 1, \dots, m $$

A cause de l'antagonisme des deux objectifs, il est impossible de trouver un minimum unique comme solution du problème. L'optimisation de la moyenne entre souvent en conflit avec la minimisation de la variance, et une décision de compromis entre les deux est nécessaire pour choisir le meilleur design. On définit alors le front de Pareto (Figure 3) qui désigne l'ensemble des solutions de conception pour lesquelles il n'existe pas de meilleure solution sur tous les objectifs. On parle aussi de solutions dominantes. La solution de conception optimale est généralement choisie comme étant le point de courbure maximale du front de Pareto. Les algorithmes génétiques, basés sur des méthodes de recherche heuristiques, sont particulièrement adaptés aux problèmes d'optimisation multiobjectif. Dans cette étude, nous utilisons l'algorithme génétique à tri non dominé de deuxième génération (NSGA-II) [6], qui est une approche évolutrice élitiste qui améliore l'exploration des solutions optimales tout en conservant une complexité raisonnable.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Illustration du front de Pareto pour un problème de conception robuste</figcaption></figure>

Ce papier explore l'extension des méthodes de couplage non-intrusives aux applications de conception robuste (Figure 4). Une première approche consiste à utiliser la méthode classique de Monte-Carlo. Cependant, l'un des principaux défis réside dans le coût des analyses stochastiques, notamment pour des modèles multiéchelles de grande taille avec un nombre important de variables aléatoires et de degrés de liberté. L'objectif est donc de développer un cadre probabiliste permettant une propagation efficace des incertitudes à travers les différentes échelles. D'une part, la méthode du chaos polynomial permet de construire des modèles de substitution capables d'estimer rapidement la réponse du système en fonction
des paramètres incertains. Plus précisément, une version non intrusive de cette méthodologie, reposant sur des techniques de régression, est appliquée [7]. Cette approche considère le modèle numérique comme une boîte noire, garantissant ainsi que l’algorithme de couplage local/global non-intrusif reste inchangé durant le processus de quantification des incertitudes. D’autre part, l’approche basée sur la méthode de perturbation permet de fournir des estimations d’incertitude de premier et de second ordre en exploitant intelligemment la structure de l’algorithme de couplage local/global. Ainsi, il devient possible de propager efficacement les incertitudes à travers plusieurs échelles et d’implémenter un processus de conception robuste en tirant parti des spécificités de l’algorithme de couplage local/global non-intrusif.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Méthodologie générale de la stratégie multiéchelle robuste envisagée</figcaption></figure>

# 4 Résultats numériques

Cette étude présente un cas pratique d'application en phase avec des problématiques réalistes rencontrées dans l'industrie aéronautique. Pour cela, nous avons élaboré un cas-test académique simple inspiré des avions servant pour le transport de charges entre divers sites de production. Notamment, nous nous intéressons au système de chargement de ce type d'avion-cargo dont la conception peut rencontrer quelques difficultés. À cause de contraintes d'assemblage, des changements de design peuvent apparaître localement et venir apporter de la souplesse supplémentaire par rapport au nominal de dimensionnement. Afin d'éviter une remise en cause tardive et coûteuse de dimensionnement, l'idée est de mettre en oeuvre une démarche de conception robuste permettant de prendre en compte dès les phases amont de conception les possibles variabilités locales. Par souci de simplicité, la structure est modélisée par un assemblage de quatre poutres en acier de mêmes caractéristiques géométriques et est discrétisée à l'aide d'éléments finis poutres de type Euler-Bernoulli. Le nouveau design est modélisé en utilisant la méthode de couplage local/global non-intrusif, où un modèle global grossier représente la conception nominale, tandis qu'un modèle local détaillé décrit la liaison flexible. La méthodologie de couplage est illustrée sur la Figure 5. L'idée est de mettre en oeuvre une démarche de dimensionnement statique à l'aide d'un modèle d'assemblage de poutres permettant de prendre en compte des variabilités de position dans les liaisons qui sont en réalité imparfaites. Ici, les poutres sont de propriétés matérielles et de longueur fixes. Ainsi, l'unique paramètre de design considéré ici est la section des poutres et l'incertitude porte sur les positions $(x_1, x_2)$ des poutres au niveau de la liaison souple qui sont modélisées par une loi uniforme.
En raison des contraintes d'assemblage, la nouvelle structure apporte une différence de flexibilité par rapport à la conception nominale. Ainsi, le problème de conception robuste consiste à déterminer la largeur $h$ de la section des poutres en acier permettant de minimiser l'écart de rigidité entre les deux modèles, malgré des variabilités de position dans la liaison souple ($x_1, x_2$) qui sont modélisées comme des paramètres stochastiques (Table 1).

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Schéma illustratif du couplage local/global non-intrusif</figcaption></figure>

<table><caption>TABLE 1 – Paramètres du problème de conception robuste</caption><tr><td>Facteurs de contrôle</td><td>$h \in [1, 100] \text{ mm}$</td></tr><tr><td>Facteurs de bruit</td><td>$x_1 \hookrightarrow U(-\frac{l}{2}, \frac{1}{2})$<br>$x_2 \hookrightarrow U(-\frac{l}{2}, \frac{1}{2})$</td></tr><tr><td>Facteurs de signal</td><td>$F = 7kN$ (dimensionnement statique)</td></tr></table>

Ici, la quantité d'intérêt est l'écart normalisé entre la flèche du modèle détaillé $U_F(h, x_1, x_2)$ et celle du modèle nominal $U_F(h_{nom}, x_1, x_2)$ au sens des moindres carrés. Le problème de conception robuste se formule comme un problème de minimisation tout en respectant les contraintes mécaniques du problème. Notamment, la contrainte équivalente de Von Mises $\sigma_{VM}$ dans la liaison souple ne doit pas dépasser la limite d'élasticité $R_e$ (avec un coefficient de sécurité $\gamma = 1.5$).

$$ \text{Minimiser :} \qquad \mu_1 [F_1(h)] , \sigma_1 [F_1(h)] \qquad (2) $$
$$ \text{Sous les contraintes : } 40 \le h \le 100 \text{ mm, } \max (\sigma_{VM}) < \frac{R_e}{\gamma} \text{ dans les boulons} $$

avec :

$$ F_1(h) = \frac{1}{2} \left( \frac{U_F(h, x_1, x_2) - U_F(h_{nom}, x_1, x_2)}{U_F(h_{nom}, x_1, x_2)} \right)^2 \qquad (3) $$

L'idée ici est de traiter le problème d'optimisation robuste précédemment défini en utilisant des méthodes alternatives de propagation des incertitudes. Ces méthodes sont plus efficaces et mieux adaptées au problème que l'approche de Monte-Carlo qui sert de solution de référence (Figure 6). La solution la plus pertinente, celle située à la courbure maximale du front de Pareto et donc représentant le meilleur compromis, peut être estimée avec moins de 5% d'erreur. De plus, contrairement à la moyenne, l'écart-type de la fonction de performance est plus difficile à estimer, particulièrement pour les petites valeurs
du paramètre de conception, c'est-à-dire lorsque le système est plus sensible aux variabilités. Il est également observé que la méthode du chaos polynomial (ou *Polynomial chaos expansion (PCE)* en anglais) fournit une approximation précise du front de Pareto tout en réduisant significativement les coûts de calcul.

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Comparaison des solutions robustes pour différentes méthodes</figcaption></figure>

Ainsi, cette méthodologie est choisie pour étudier l'impact du couplage des algorithmes de couplage local/global non-intrusifs sur l'optimisation robuste. L'influence de la convergence de l'algorithme local/global non-intrusif sur l'estimation du front de Pareto est analysée (Figure 7). Dès que l'erreur en résidu de l'algorithme de couplage local/global non-intrusif est suffisamment faible, les fronts de Pareto estimés convergent vers celui de référence et restent admissibles. Bien que les solutions robustes ne soient pas entièrement optimales, elles permettent de réaliser une phase de pré-conception robuste et rapide. Aussi, il convient de noter que la méthode de *submodeling* ne donne pas de résultats satisfaisants, en raison principalement de la nature localisée des incertitudes.

# 5 Conclusion

En conclusion, le développement d'une approche de conception robuste pour la simulation multiéchelle représente un pas important vers la gestion efficace des incertitudes en ingénierie. Pour cela, les performances de différentes méthodes de propagation d'incertitudes (chaos polynomial, méthode par perturbation) ont été analysées et comparées sur un cas-test académique afin d'aboutir à une approche efficace en fonction du cadre d'étude. Par la suite, afin de gérer efficacement les incertitudes, il est intéressant d'associer à l'approche multiéchelle une stratégie de modélisation hybride dans lequel le modèle global est couplé à un métamodèle issu d'une modélisation haute-fidélité et stochastique, permettant ainsi des réponses en temps court dans un cadre de calcul multi-requêtes. Pour mener l'exploration rapide du domaine de conception potentiel et obtenir des réponses rapides dans le test de diverses configurations, on envisage d'employer localement les méthodes de réduction de modèle. Ces approches permettront de trouver un équilibre entre précision et coût de calcul, renforçant ainsi l'efficacité des méthodologies de conception robuste pour les problèmes complexes en ingénierie.
<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Comparaison des fronts de Pareto pour différentes convergences de l’algorithme local/global</figcaption></figure>

Références

[1] C. Zang, M. I. Friswell, and J. E. Mottershead. A review of robust optimal design and its application in dynamics. *Computers & Structures*, 83(4-5):315-326, 2005.

[2] M. Guedri, S. Cogan, and N. Bouhaddi. Robustness of structural reliability analyses to epistemic uncertainties. *Mechanical Systems and Signal Processing*, 28:458-469, 2012.

[3] M. Ouisse and S. Cogan. A Decision-making Methodology for the Robust Design of Spot Welds in Automotive Structures. *IMAC XXVI*, 2008.

[4] L. Gendre, O. Allix, P. Gosselet, and F. Comte. Non-intrusive and exact global/local techniques for structural problems with local plasticity. *Computational Mechanics*, 44:233-245, 2009.

[5] P. Minigher, A. Arteiro, A. Turon, J. Fatemi, S. Guinard, L. Barrière, and P.P. Camanho. On an efficient global/local stochastic methodology for accurate stress analysis, failure prediction and damage tolerance of laminated composites. *International Journal of Solids and Structures*, 303:113026, 2024.

[6] Kalyanmoy Deb. A Fast and Elitist Multiobjective Genetic Algorithm : NSGA-II. *IEEE transactions on evolutionary computation*, 6(2):182-197, 2002.

[7] M. Berveiller, B. Sudret, and M. Lemaire. Stochastic finite element : a non intrusive approach by regression. *European Journal of Computational Mechanics/Revue Européenne de Mécanique Numérique*, 15(1-3):81-92, 2006.