La procédure TOPOPTIM (v3.0)

La procédure TOPOPIM de Cast3M reprend les algotithmes 99 et 88 lignes de [SIGMUND-2001] et [ANDREASSEN-AL-2011] :

  • Le problème d'optimisation est celui de la minimisation de la compliance sous contrainte d'une fraction volumique imposée

  • La méthode de pénalisation SIMP est utilisée

  • L'optimisation est réalisée par critère d'optimalité et le multiplicateur de Lagrange (associé à la contrainte de fraction volumique) est calculé par dichotomie

  • Le filtrage de la sensibilité est réalisé par convolution

  • L'utilisateur peut imposer des zones figées qui feront obligatoirement partie de la topologie optimisée \((x_e=1)\)

De surcroît, elle élargie son cadre d'application à différents types de modélisation (2D/3D, maillages linéaires/quadratiques, tout type d'éléments, etc...) et dispose de fonctionnalités supplémentaires :

  • L'utilisation de maillages non structurés et contenant plusieurs types d'éléments

  • L'application à des comportements mécaniques non linéaires (plasticité, contact, grands déplacements ...)

  • La prise en compte de plusieurs cas de chargements. L'utilisateur peut définir un nombre \(N_c\) de problèmes à résoudre et minimiser une fonctionnelle sommant les compliances

  • L'application à des problèmes thermiques et/ou thermo mécaniques. Dans ce second cas, les compliances mécaniques \(f^{\textrm{méca}}\) et thermiques \(f^{\textrm{ther}}\) sont pondérées et additionnées

  • La possibilité d'utiliser des restrictions géométriques avec des symétries (centrale, axiale et plane) ou des conditions de pérodicité (axiale ou circulaire)

  • La possibilité de faire évoluer les paramètres au cours des itérations de l'optimisation

  • La prise en compte de synthèse de mécanisme souple

Quelques définitions

La fonction objectif pour un problème mécanique est la compliance qui correspond au travail des forces extérieures, ou aussi à l'énergie de déformation :

(1)\[f^{\textrm{méca}}(\textbf{x}) = \textbf{U}^T(\textbf{x}).\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) = \int_{\Omega} \sigma(\textbf{x}):\varepsilon(\textbf{x}) dV\]

La fonction objectif pour un problème thermique est construite de manière similaire à celle de la mécanique :

(2)\[f^{\textrm{ther}}(\textbf{x}) = \Theta^T(\textbf{x}).\mathbfcal{\Lambda}(\textbf{x}).\Theta(\textbf{x}) = \int_{\Omega} \phi(\textbf{x}).\textrm{grad}\,\Theta(\textbf{x}) dV\]

avec :

  • \(\Theta\) le vecteur des températures aux noeuds du maillage

  • \(\mathbfcal{\Lambda}\) la matrice de conductivité globale, assemblée sur le maillage

  • \(\lambda\) la conductivité thermique

  • \(\textrm{grad}\, \Theta\) le gradient de température dans les éléments

  • \(\phi = -\lambda\,\textrm{grad}\, \Theta\) la densité de flux thermique dans les éléments

Étapes de la procédure

La procédure TOPOPIM est segmentée en plusieurs étapes avec des sous procédures spécifiques. Outre la facilité de la compréhension et de la maintenance, cette segmentation permet à l'utilisateur/développeur de pouvoir modifier la procédure facilment pour l'adapter à son propre problème d'optimisation.

Les étapes sont les suivantes :

  • Pré traitement (procédure TOPOBOOT)

    On initialise une table de travail (indice WTABLE) contenant les indicateurs et les valeurs des paramètres choisis par l'utilisateur ainsi que des variables de travail (champs unitaires, volumes élémentaires, etc). On initialise notament la topologie \(\textbf{x}\) comme uniforme et égale à la fraction volumique imposée \(x_e=f\)

  • Informations (procédure TOPOINFO)

    Il s'agit d'afficher les informations principales du calcul.

  • Début de la boucle d'optimisation :

    1. Évolution des paramètres (procédure TOPOFCTR)

      On fait évoluer les principaux coefficients (d'amortissement, de pénalisation et de niveau de gris) au cours des itérations, à condition que l'utilisateur ai fourni une liste de valeurs pour ces coefficients en fonction du numéro d'itération.

    2. Pénalisation des densités (procédure TOPODENS)

      On calcule les champs de densité pénalisées \(\tilde{\textbf{x}}\) selon la topologie courante \(\textbf{x}\), le coefficient de pénalisation \(p\) et la raideur minimale \(E_{\textrm{min}}\) et/ou la conductivité minimale \(\lambda_{\textrm{min}}\) :

      (3)\[\begin{split}\tilde{x}_e^{\textrm{méca}}& =\frac{E_{\textrm{min}}}{E_0}+\left(1-\frac{E_{\textrm{min}}}{E_0}\right)x_e^p \\ \tilde{x}_e^{\textrm{ther}}& =\frac{\lambda_{\textrm{min}}}{\lambda_0}+\left(1-\frac{\lambda_{\textrm{min}}}{\lambda_0}\right)x_e^p\end{split}\]
    3. Mise à jour des zones actives et matériau (procédure TOPOACTI)

      On met à jour les zones actives du maillage et des modèles. Il s'agit des éléments où la densité pénalisée \(\tilde{x}_e\) est supérieure à un seuil d'activation \(x_a\). Les modèles et caractéristiques sont alors réduits sur le maillage actif. Les caractéristiques matériau sont également mises à jour selon la densité pénalisée avec la procédure TOPOMATE :

      (4)\[\begin{split}\tilde{E}_e & = E_0 \tilde{x}_e^{\textrm{méca}} \\ \tilde{\lambda}_e & = \lambda_0 \tilde{x}_e^{\textrm{ther}}\end{split}\]
    4. Résolution du problème physique (procédure TOPORESO).

      La résolution du problème mécanique et/ou thermique se fait soit avec le solveur RESO (problème linéaire), soit avec PASAPAS (problème non linéaire). S'il y a plusieurs cas de chargements, les \(N_c\) problèmes sont tous résolus successivement.

    5. Instructions personnelles (optionnel, procédure TOPOPERS)

      Il s'agit d'un point de branchement donné à l'utilisateur pour faire appel à une procédure personnelle et ajouter des instructions supplémentaires, ou bien remplacer une procédure et ainsi modifier l'algorithme.

    6. Objectif et sensibilités (procédure TOPOSENS)

      On calcule la fonction objectif \(f(\tilde{\textbf{x}})\) ainsi que le champ de sensibilité \(\frac{\partial f}{\partial \tilde{x}_e}\). S'il y a un plusieurs cas de chargement, les \(N_c\) objectifs/sensibilités sont sommés. Dans le cas thermo mécanique, les complicances des deux physiques sont donc aussi sommées mais peuvent être pondérées différement par les coefficients \(\omega^{\textrm{méca}}\) et \(\omega^{\textrm{ther}}\) :

      (5)\[f(\tilde{\textbf{x}}) = \frac{1}{N_c} \sum_{i=1}^{N_c} \left( \omega^{\textrm{méca}} f_i^{\textrm{méca}}(\tilde{\textbf{x}}) + \omega^{\textrm{ther}} f_i^{\textrm{ther}}(\tilde{\textbf{x}}) \right)\]
      (6)\[\frac{\partial f}{\partial \tilde{x}_e} = \frac{1}{N_c} \sum_{i=1}^{N_c} \left( \omega^{\textrm{méca}} \frac{\partial f_i^{\textrm{méca}}}{\partial \tilde{x}_e} + \omega^{\textrm{ther}} \frac{\partial f_i^{\textrm{ther}}}{\partial \tilde{x}_e} \right)\]
    7. Restrictions géométriques (procédure TOPORSTR)

    8. Filtrage (procédure TOPOFILT)

      Le filtrage de la sensibilité peut être réalisé selon 2 méthodes :

      • Le filtre GIBIANE qui procède par des applications successives de l'opérateur CHAN 'CHPO' puis CHAN 'CHAM' sur le champ \(\tilde{x}_e \frac{\partial f}{\partial \tilde{x}_e}\) ayant pour effet de le lisser

      • Le filtre MATRICE qui applique une convolution (décrite ici) sur le champ \(\tilde{x}_e \frac{\partial f}{\partial \tilde{x}_e}\) via l'opérateur MFIL et correspond au filtrage habituellement utilisé, notamment dans [SIGMUND-2001]

    9. Optimisation de la topologie (procédure TOPOLOGY)

      La mise à jour de la topologie \(\textbf{x}\) s'effectue selon un algorithme de critère d'optimalité (décrit ici) où le calcul du multiplicateur de Lagrange est fait par dichotomie (décrit ici).

    10. Sauvegarde des résultats (procédure TOPOSAUV)

    11. Tracé et infos éventuel de la topologie et affichage des informations sur l'itération courante

    12. Test de convergence

      La boucle d'optimisation est quittée si l'incrément maximal de densité entre deux itérations est inférieur au critère \(Z_{\textrm{stop}}\), ou bien si le nombre maximal d'itérations \(N_{\textrm{it}}\) est atteint.

  • Fin de la boucle d'optimisation

Variables/paramètres et indices de la table de calcul

Le tableau ci-dessous fait la correspondance entre les variables/paramètres des problèmes d'optimisation présentées ici et les indices de la table de calcul utilisée par TOPOPTIM. Certaines grandeurs sont stockées directement dans la table car elles sont fixes ou bien sont des résultats de calcul rendus à l'utilisateur, d'autres sont stockées dans la sous table WTABLE car elles sont mises à jour à chaque itération et sont des intermédiaires de calcul temporaires.

Les valeurs choisies par défaut de certains paramètres sont également indiquées. Ce tableau n'est pas une liste exhaustive de la table, le lecteur intéressé peut consulter le code de la procédure TOPOBOOT qui initialise cette table.

Variables/paramètres et indices de la table de calcul

Variable / Paramètre + Rôle

Indice dans la table

Valeur par défaut

\(x_e\)

Variable de conception

WTABLE . TOPOLOGIE

\(\tilde{x}_e^{\textrm{méca}}\)

Densité (pour la mécanique)

WTABLE . MECANIQUE . DENSITE

\(\tilde{x}_e^{\textrm{ther}}\)

Densité (pour la thermique)

WTABLE . THERMIQUE . DENSITE

\(f\)

Fraction volumique cible

FRACTION_VOLUME

0,4

\(p\)

Coefficient de pénalisation

WTABLE . FACTEUR_P

3

\(\eta\)

Coefficient d'amortissement

WTABLE . FACTEUR_D

0,5

\(x_\textrm{min}\)

Densité minimale

TOPOLOGIE_MIN

0

\(m\)

Limite d'incrément de densité

TOPOLOGIE_MAX_INC

0,2

\(\frac{E_{\textrm{min}}}{E_0}\)

Ratio de raideur mécanique

RAPPORT_RAIDEURS_MECANIQUES

10-8

\(\frac{\lambda_{\textrm{min}}}{\lambda_0}\)

Ratio de raideur thermique

RAPPORT_RAIDEURS_THERMIQUES

10-3

\(\omega^{\textrm{méca}}\)

Poids de mécanique

POIDS_ENERGIE_DEFO

1

\(\omega^{\textrm{ther}}\)

Poids de la thermique

POIDS_TEMPERATURE

1

\(N_c\)

Nombre de cas de chargement

WTABLE . NB_CAS

1

\(x_a\)

Seuil d'activation des éléments

SEUIL

10-9

\(r_\textrm{min}\)

Rayon de filtrage

FILTRE_RAYON

\(q\)

Exposant pour le filtrage

FILTRE_EXPOSANT

1

\(f\)

Fonction objectif

WTABLE . OBJECTIF

\(\dfrac{\partial f}{\partial \tilde{x}_e}\)

Sensibilité de la fonction objectif

WTABLE . SENSIBILITE

\(Z_{\textrm{stop}}\)

Critère de convergence

CRITERE

0,01

\(N_{\textrm{it}}\)

Nombre max. d'itérations

MAX_CYCLES

100