.. _sec:opti_topo_appli: Illustration sur un cas mécanique ================================= Une mise en donnée de :ref:`l'algorithme d'optimisation précédent ` est fournie en :ref:`annexe ` et disponible sur le site Cast3M dans `top_oc.dgibi `_. On propose ci-dessous une analyse détaillée des variables calculées et leur lien avec les éléments théoriques précedemment énoncés. Cet exemple Gibiane reproduit le cas test de la *poutre MBB* des articles de référence [SIGMUND-2001]_ et [ANDREASSEN-AL-2011]_. Définition du problème d'optimisation et paramètres --------------------------------------------------- Il s'agit d'optimiser la topologie d'une **poutre en flexion 3 points**. Le problème possède des symmétries et est modélisé en hypohèse 2D contraintes planes, en ne représentant qu'une moitié de la poutre. Celle-ci est soumise à une force ponctuelle verticale en haut à gauche, supportée verticalement en bas à droite et soumise à des conditions de symmétrie sur le bord gauche. .. figure:: figures/99_lines_img_1.png :width: 15cm :align: center Pour l'optimisation, on choisit : - une **fonction objectif** : la compliance :math:`f(\textbf{x}) = \textbf{U}^T.\textbf{F}` - une **fonction contrainte** sur la fraction volumique : :math:`g(\textbf{x}) = \phi(\textbf{x}) - \phi_0 = 0` avec :math:`\phi_0 = 50\%` - les dimensions longueur :math:`l=60` m, hauteur :math:`h=20` m et épaisseur :math:`e=1` m - un module d'Young unitaire : :math:`E_0=1` Pa, :math:`\nu=0,3` - un chargement de force unitaire - un rayon de filtr : :math:`r_{\textrm{min}}=0,04.l` - des paramètres d'optimisation : :math:`p=3`, :math:`\eta=0,5`, :math:`m=0,2` et :math:`E_{\textrm{min}}=10^{-9}.E_0` Dans le programme, la variable ``ft`` permet de choisir le type de filtre appliqué : - ``ft = 1`` pour utiliser un filtre sur la sensibilité :math:`\frac{\partial f}{\partial x_e}` - ``ft = 2`` pour utiliser un filtre sur les densités - ``ft = 3`` pour utiliser un filtre sur les densité et un seuillage par fonction de Heaviside Initialisations --------------- On commence par créer qelques objets utiles pour la suite : - un champ unitaire par élément : ``un`` - un champ matériau unitaire : ``maun`` - le champ des volumes élémentaires ``ve``, obtenu grâce à l'opérateur `INTG 'ELEM' `_ en intégrant le champ unitaire ``un`` - le volume du domaine de conception : ``vtot`` - une mesure de l'incrément de topologie : ``change``, qui servira de critère d'arrêt Puis on initialise la topologie ``x`` comme un champ homogène :math:`x_e=\phi_0` afin de satisfaire la contrainte de volume. .. admonition:: Initialisation : champs utiles et topologie initiale .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 104-112 Les champs des autres densités, filtrés, seuillées, sont calculées selon l'ioption ``ft`` choisie. .. admonition:: Initialisation : densités filtrées et seuillées initiales .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 114-121 On calcule la matrice ``kfil``, réalisant le :ref:`filtrage ` grâce aux opérateurs `MFIL `_ (pour le numérateur) et `NFIL `_ (pour normer le résultat par le dénominateur). Notons que pour cela, il est nécessaire de disposer du maillage ``ptg`` des centres de gravité du maillage ainsi que du champ par points ``wg`` des volumes :math:`V_e` de chaque éléments, exprimé sur ces centres de gravité. .. admonition:: Initialisation : matrice de filtrage .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 126-133 Pénalisation et résolution du problème mécanique ------------------------------------------------ On démarre ensuite une boucle d'optimisation, limitée à 500 itérations. On calcule le champ ``ep``, correspondant au champ de module d'Young pénalisé :math:`E(x_e)`, puis ``k`` la matrice de rigidité pénalisée :math:`\mathbfcal{K}(\textbf{x})`, à laquelle on adjoint la matrice des blocages ``blo`` pour former la rigidité totale ``kbc``. On résoud ensuite le problème mécanique :math:`\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) =\textbf{F}` en calculant les déplacements ``u`` avec l'opérateur `RESO `_. .. admonition:: Pénalisation de la rigidité et résolution .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 145,147-153 Calcul des fonctions objectif, contrainte et de leurs sensibilités ------------------------------------------------------------------ On peut évaluer ``c``, la fonction objectif :math:`f(\textbf{x})`, via l'opérateur de produit scalaire `PSCA `_. Le champ ``dc`` de sensibilité :math:`\frac{\partial f}{\partial x_e}` est calculé en passant par les champs de déformation ``eps`` et de contraintes ``sigun``. En effet, si la compliance est égale au travail des forces extérieures, elle est aussi égale au travail des efforts intérieurs. La sensibilité de :math:`f` peut s'écrire : .. math:: :name: eq:opti_topo_sensibilite_3 \frac{\partial f}{\partial x_e} &= -p(x_e)^{p-1}(E_0-E_{\textrm{min}}) \textbf{u}_e^T.\mathbfcal{k}_0.\textbf{u}_e\\ &= -p(x_e)^{p-1}(E_0-E_{\textrm{min}}) \int_{V_e} \pmb{\varepsilon}:\mathbfcal{C}_0:\pmb{\varepsilon}~dV où :math:`\pmb{\varepsilon}` est le tenseur de déformation ``eps`` et :math:`\mathbfcal{C}_0` la matrice de Hooke pour un module d'Young unitaire. Le champ d'énergie ``eneun`` (le double produit contracté :math:`\pmb{\varepsilon}:\mathbfcal{C}_0:\pmb{\varepsilon}`) est obtenu grâce à l'opérateur `ENER `_ et son intégrale par élément avec `INTG 'ELEM' `_. .. admonition:: Calcul de la fonction objectif et de sa sensibilité .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 154-162 La fonction contrainte ``g`` est évaluée simplement à partir de la fraction volumique ``vfx``, tout comme le champ ``dg`` de sensibilité :math:`\frac{\partial g}{\partial x_e}=\frac{V_e}{V_0}` : .. admonition:: Calcul de la fonction contrainte et de sa sensibilité .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 163-167 Procédure de filtrage --------------------- :ref:`L'opération de filtrage ` est codée de manière générique dans un procédure ``HFILT`` en début de programme, pouvant s'appliquer à n'importe quel champ ``cham1``. Elle consiste à transformer ce champ en un vecteur ``chp1``, le multiplier par la matrice de filtrage ``kfil`` pour obtenir le vecteur filtré ``chp2`` et reconstruie un champ filtré ``cham2`` sur le même support (centres de gravité des éléments). .. admonition:: Filtrage de la sensibilité .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 37-42 Filtrage de la sensibilité (``ft=1``) ------------------------------------- Si l'on :ref:`filtre la sensibilité `, on applique la procédure de filtrage ``HFILT`` au produit de champ ``x*dc``, que l'on divise ensuite par :math:`\textrm{max}(10^{-3},x_e)`, gr^ce à l'opérateur `BORN `_. .. admonition:: Filtrage de la sensibilité .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 168-171 Filtrage des densités (``ft=2`` ou ``ft=3``) -------------------------------------------- *A faire* Seuillage des densités (``ft=3``) --------------------------------- *A faire* Optimisation par critère d'optimalité ------------------------------------- La mise à jour de la topologie (passage du champ ``x`` à ``xnew``) suivant le :ref:`schéma ` présenté précédemment est réalisée en suivant :ref:`l'algorithme de dichotomie ` pour la recherche du multiplicateur de Lagrange ``lmid`` qui nécessite une nouvelle boucle (limitée à 200 itérations). La limitation d'incrément :math:`m` et le recpect des bornes :math:`0 \le x_e \le 1` sont réalisées grâce aux opérateurs `BORN `_ et `MASQ `_. .. admonition:: Optimisation (critère d'optimalité) .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 191-206 *A faire : filtrage des densités* La vérification de la contrainte est faite en calculant ``vfx``, la fraction volumique :math:`\phi(\textbf{x})`, et la nouvelle valeur ``g`` de la fonction contrainte :math:`g(\textbf{x})=0`. La division de l'intervalle :math:`[\mathcal{L}^-,\mathcal{L}^+]` est faite selon le signe de ``g`` .. admonition:: Optimisation (critère d'optimalité) .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 218-226 Un bilan de l'itération est fait, puis un cirtère d'arrêt de la boucle d'optimisation est proposé lorsque l'incrément maximal de densité ``change`` devient inférieur à ``changmax`` (choisit égal à 0,01). .. admonition:: Fin de boucle et critère d'arrêt .. literalinclude:: dgibi/top_oc.dgibi :language: gibiane :lines: 227-228,233-234,243-244,247-249 Les résultats de cette optimisation sont présentés dans l'animation ci-dessous qui montre les topologies (champs par éléments de densités) obtenues au cours des itérations pour ``ft = 2`` (filtrage des densités sans seuillage). .. figure:: figures/top_oc.gif :name: fig:top_anim1 :width: 20cm :align: center Animation des topologies au cours de l'optimisation (déformée x 1000)