Illustration sur un cas mécanique

Une mise en donnée de l'algorithme d'optimisation précédent est fournie en 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.

../_images/99_lines_img_1.png

Pour l'optimisation, on choisit :

  • une fonction objectif : la compliance \(f(\textbf{x}) = \textbf{U}^T.\textbf{F}\)

  • une fonction contrainte sur la fraction volumique : \(g(\textbf{x}) = \phi(\textbf{x}) - \phi_0 = 0\) avec \(\phi_0 = 50\%\)

  • les dimensions longueur \(l=60\) m, hauteur \(h=20\) m et épaisseur \(e=1\) m

  • un module d'Young unitaire : \(E_0=1\) Pa, \(\nu=0,3\)

  • un chargement de force unitaire

  • un rayon de filtr : \(r_{\textrm{min}}=0,04.l\)

  • des paramètres d'optimisation : \(p=3\), \(\eta=0,5\), \(m=0,2\) et \(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é \(\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 \(x_e=\phi_0\) afin de satisfaire la contrainte de volume.

Initialisation : champs utiles et topologie initiale

** Volume of each element (ve) and total volume (vtot)
un       = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
ve       = INTG 'ELEM' mod un maun ;
vtot     = INTG mod un maun ;

** Initial density field
xini     = volfrac ;
x        = MANU 'CHML' mod 'SCAL' xini 'GRAVITE' ;
change   = 1. ;

Les champs des autres densités, filtrés, seuillées, sont calculées selon l'ioption ft choisie.

Initialisation : densités filtrées et seuillées initiales

** Physical density field
SI ((ft EGA 1) OU (ft EGA 2)) ;
  xphys    = x ;
FINSI ;
SI (ft EGA 3) ;
  xtilde   = x ;
  xphys    = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
FINSI ;

On calcule la matrice kfil, réalisant le 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 \(V_e\) de chaque éléments, exprimé sur ces centres de gravité.

Initialisation : matrice de filtrage

** Gravity centers and filtering matrix
*  the weight depends on the volume of each element
ptg      = un POIN 'SUPERIEUR' 0. ;
wg       = MANU 'CHPO' ptg 'SCAL' (EXTR ve 'VALE' 'SCAL') ;
kfil     = MFIL wg rmin 1. 0. ;
ung      = MANU 'CHPO' ptg 1 'SCAL' 1. ;
ks       = kfil * ung ;
kfil     = NFIL kfil ks ;

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é \(E(x_e)\), puis k la matrice de rigidité pénalisée \(\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 \(\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x}) =\textbf{F}\) en calculant les déplacements u avec l'opérateur RESO.

Pénalisation de la rigidité et résolution

REPE b1 500 ;
* penalization of the stiffness matrix (modified SIMP)
  ep       = emin + ((xphys ** penal) * (e0 - emin)) ;
  map      = MATE mod 'YOUN' ep 'NU' nu 'DIM3' e ;
  k        = RIGI mod map ;
* resolution of the FE problem
  kbc      = k ET blo ;
  u        = RESO kbc f ;

Calcul des fonctions objectif, contrainte et de leurs sensibilités

On peut évaluer c, la fonction objectif \(f(\textbf{x})\), via l'opérateur de produit scalaire PSCA.

Le champ dc de sensibilité \(\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 \(f\) peut s'écrire :

(1)\[\begin{split}\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\end{split}\]

où \(\pmb{\varepsilon}\) est le tenseur de déformation eps et \(\mathbfcal{C}_0\) la matrice de Hooke pour un module d'Young unitaire. Le champ d'énergie eneun (le double produit contracté \(\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'.

Calcul de la fonction objectif et de sa sensibilité

* objective function: compliance = u^T.F
* value
  c        = MAXI (RESU (PSCA u f (MOTS 'UX' 'UY') (MOTS 'FX' 'FY'))) ;
* sensitivity (gradient with respect to the physical variable)
  eps      = EPSI 'LINE' mod u ;
  sigun    = ELAS mod maun eps ;
  eneun    = ENER mod eps sigun ;
  eneun    = INTG 'ELEM' mod eneun ;
  dc       = -1. * penal * (xphys ** (penal - 1.)) * (e0 - emin) * eneun ;

La fonction contrainte g est évaluée simplement à partir de la fraction volumique vfx, tout comme le champ dg de sensibilité \(\frac{\partial g}{\partial x_e}=\frac{V_e}{V_0}\) :

Calcul de la fonction contrainte et de sa sensibilité

* constraint function: g = vf(x) - volfrac
* value
  g        = vfx - volfrac ;
* sensitivity (gradient with respect to the physical variable)
  dg       = ve / vtot ;

Procédure de filtrage

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).

Filtrage de la sensibilité

** Filtering procedure
DEBP HFILT cham1*'MCHAML' mo*'MMODEL' mf*'RIGIDITE' mpt*'MAILLAGE' ;
  chp1     = MANU 'CHPO' mpt 'SCAL' (EXTR cham1 'VALE' 'SCAL') ;
  chp2     = mf * chp1 ;
  cham2    = MANU 'CHML' mo 'REPA' 'SCAL' (EXTR chp2 'VALE' mpt) 'TYPE' 'SCALAIRE' 'GRAVITE' ;
FINP cham2 ;

Filtrage de la sensibilité (ft=1)

Si l'on filtre la sensibilité, on applique la procédure de filtrage HFILT au produit de champ x*dc, que l'on divise ensuite par \(\textrm{max}(10^{-3},x_e)\), gr^ce à l'opérateur BORN.

Filtrage de la sensibilité

* ft = 1  --> filtering the complicance sensitivity
  SI (ft EGA 1) ;
    dc       = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
  FINSI ;

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 schéma présenté précédemment est réalisée en suivant 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 \(m\) et le recpect des bornes \(0 \le x_e \le 1\) sont réalisées grâce aux opérateurs BORN et MASQ.

Optimisation (critère d'optimalité)

* update the topology by optimization
  l1       = 0. ;
  l2       = 1.E9 ;
  REPE b2 200 ;
    SI (((l2 - l1) / (l1 + l2)) < 0.001) ;
      QUIT b2 ;
    FINSI ;
    lmid     = 0.5 * (l1 + l2) ;
    b        = -1. * dc / (lmid * dg) ;
    xinf     = BORN (x - move) 'MINIMUM' xmin ;
    xsup     = BORN (x + move) 'MAXIMUM' xmax ;
    xnew     = x * (b ** 0.5) ;
    minf     =  (xnew - xinf) MASQ 'INFERIEUR' 0. ;
    mmil     = ((xnew - xinf) MASQ 'SUPERIEUR' 0.) * ((xnew - xsup) MASQ 'INFERIEUR' 0.) ;
    msup     =  (xnew - xsup) MASQ 'SUPERIEUR' 0. ;
    xnew     = (xinf * minf) + (xnew * mmil) + (xsup * msup) ;

A faire : filtrage des densités

La vérification de la contrainte est faite en calculant vfx, la fraction volumique \(\phi(\textbf{x})\), et la nouvelle valeur g de la fonction contrainte \(g(\textbf{x})=0\). La division de l'intervalle \([\mathcal{L}^-,\mathcal{L}^+]\) est faite selon le signe de g

Optimisation (critère d'optimalité)

*   updating the volume fraction and the constraint function
    vfx      = (INTG mod xphys maun) / vtot ;
    g        = vfx - volfrac ;
    SI (g > 0.) ;
      l1       = lmid ;
    SINON ;
      l2       = lmid ;
    FINSI ;
  FIN b2 ;

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).

Fin de boucle et critère d'arrêt

* summary of the current iteration
  change   = MAXI 'ABS' (xnew - x) ;
* preparing the next iteration
  x        = xnew ;
* stop criterion
  SI (change < changmax) ;
    QUIT b1 ;
  FINSI ;
FIN  b1 ;

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).

../_images/top_oc.gif

Animation des topologies au cours de l'optimisation (déformée x 1000)