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

Définition du problème d'optimisation

Il s'agit d'optimiser la poutre en flexion présentée plus haut

../_images/99_lines_img_1.png

Pour l'optimisation, on choisit :

  • une fonction objectif : la compliance \(\psi(\textbf{x}) = \textbf{U}^T(\textbf{x}).\mathbfcal{K}(\textbf{x}).\textbf{U}(\textbf{x})\)

  • une contrainte sur le volume : \(f=40\%\) du domaine de conception

  • les paramètres d'optimisation : \(p=3\), \(\eta=0,5\), \(m=0,1\) et \(x_{\textrm{min}}=0,001\).

Initialisations

On initialise la topologie x avec des densités homogènes \(x_e=f\) afin de satisfaire la contrainte de volume. Le volume cible est nommé vcib.

On calcule la matrice de filtrage kfil, intervenant dans l'équation (8) avec l'opérateur MFIL. Notons que pour cela, il est nécessaire de disposer du maillage mcg 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é. Le champ des volumes élémentaires vole est obtenu grâce à l'opérateur INTG 'ELEM' en intégrant le champ un unitaire par élément.

Les volumes utiles sont aussi calculés : v0 le volume du domaine de conception, vx le volume de la topologie x courante et vcib le volume cible.

Initialisation : topologie initiale et matrice de filtrage

46** Initialisation de la topologie (avec la fraction volumique cible)
47x      = MANU 'CHML' mod 'SCAL' fv 'GRAVITE' ;
48
49** Matrice de filtrage
50un     = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
51vole   = INTG mod un 'ELEM' ;
52mcg    = un POIN 'SUPERIEUR' 0. ;
53wg     = PROI mcg vole ;
54kfil   = MFIL wg rmin 1. 0. ;
55
56** Volume plein, initial et cible
57v0     = INTG mod un ;
58vx     = INTG mod x ;
59fvx    = vx / v0 ;
60vcib   = v0 * fv ;
61chgx   = 0. ;

On démarre ensuite une boucle d'optimisation limitée à 100 itérations.

Pénalisation et résolution du problème mécanique

On calcule alors rip, la matrice de rigidité pénalisée \(\mathbfcal{K}(\textbf{x})\) de la topologie courante selon la loi puissance de la méthode SIMP. Le comportement étant isotrope, le module d'Young pénalisé yop de chaque élément vaut \(E_e=(x_e)^pE_0\) avec \(E_0\) le module d'Young du matériau.

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

63** Boucle d'optimisation topologique
64liso   = PROG 0. 'PAS' 0.05 1. ;
65REPE b1 100 ;
66* pénalisation de la rigidité
67  yop    = (x ** p) * yo ;
68  map    = MATE mod 'YOUN' yop 'NU' nu ;
69* résolution du problème mécanique
70  rip    = RIGI mod map ;
71  u      = RESO (rip ET blo) f ;

Calcul de la fonction objectif et des sensibilités

On peut calculer la valeur psi de la fonction objectif \(\psi(\textbf{x}) = \textbf{F}^T(\textbf{x}).\textbf{U}(\textbf{x})\) en remarquant que si celle-ci est égale au travail des forces extérieures, elle est donc aussi égale au travail des efforts intérieurs et peut donc s'obtenir par :

\[\psi(\textbf{x}) = \int_{\Omega} \sigma(\textbf{x}):\varepsilon(\textbf{x}) dV\]

où \(\sigma\) et \(\varepsilon\) désignent les champs de contraintes et déformations sig et eps. Le champ du double produit contracté \(\sigma:\varepsilon\) est obtenu grâce à l'opérateur ENER et son intégrale par INTG.

Le champ dpsi de sensibilité de la fonction objectif (6) s'exprime alors en fonction de la matrice de Hooke \(\mathbfcal{C}_0\) du matériau plein :

\[\frac{\partial \psi}{\partial x_e} = -p(x_e)^{p-1} \varepsilon^T(x_e).\mathbfcal{C}_0.\varepsilon(x_e)\]

Calcul de la fonction objectif et de sa sensibilités

72* fonction objectif : compliance = uT.K.u = Int(sig:eps)dV
73  eps    = EPSI 'LINE' mod u ;
74  sig    = ELAS mod map eps ;
75  psi    = INTG mod (ENER mod sig eps) ;
76* sensibilité
77  sig0   = ELAS mod ma0 eps ;
78  ene0   = ENER mod eps sig0 ;
79  ene    = CHAN ene0 mod 'GRAVITE' ;
80  dpsi   = (-1. * p * (x ** (p - 1.)) * ene) ;

Filtrage

L'étape de filtrage de la sensibilité est réalisée en multipliant la matrice de filtrage kfil par le champ par point xdpsi = x * dpsi représentant le produit \(x_f\dfrac{\partial \psi}{\partial x_f}\) dans l'équation (8).

Filtrage de la sensibilité

81* filtrage de la sensibilité
82  xdpsi  = x * dpsi ;
83  xdpsi0 = PROI mcg xdpsi ;
84  xdpsi1 = kfil * xdpsi0 ;
85  xdpsi  = MANU 'CHML' mod 'REPA' 'SCAL' (EXTR xdpsi1 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
86  dpsi   = xdpsi / x ;

Optimisation par critère d'optimalité

La mise à jour de la topologie (passage du champ x à xnew) suivant le schéma (4) 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 à 100 itérations).

La limitation d'incrément \(m\) et le recpect des bornes \(x_\textrm{min} \le x_e \le 1\) sont réalisées grâce aux opérateurs BORN et MASQ.

La vérification de la contrainte de volume est faite en calculant le volume vxnew de chaque toplogie xnew et en le comparant au volume cible vcib.

Optimisation (critère d'optimalité)

 94* optimisation d'une nouvelle topologie (méthode du critère d'optimalité)
 95  l1     = 0. ;
 96  l2     = 100000000. ;
 97  REPE b2 100 ;
 98    SI ((l2 - l1) < 0.0001) ;
 99      QUIT b2 ;
100    FINSI ;
101    lmid   = (l1 + l2) / 2. ;
102    b      = -1. * dpsi / (lmid * vole) ;
103    xinf   = BORN (x - m) 'MINIMUM' xmin ;
104    xsup   = BORN (x + m) 'MAXIMUM' xmax ;
105    xnew   = x * (b ** eta) ;
106    minf   =  (xnew - xinf) MASQ 'INFERIEUR' 0. ;
107    mmil   = ((xnew - xinf) MASQ 'SUPERIEUR' 0.) * ((xnew - xsup) MASQ 'INFERIEUR' 0.) ;
108    msup   =  (xnew - xsup) MASQ 'SUPERIEUR' 0. ;
109    xnew   = (xinf * minf) + (xnew * mmil) + (xsup * msup) ;
110    vxnew  = INTG mod xnew ;
111    SI (vxnew > vcib) ;
112      l1     = lmid ;
113    SINON ;
114      l2     = lmid ;
115    FINSI ;
116  FIN b2 ;

Un affichage 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é chgx est inférieur à 0,01

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

117* bilan de l'itération
118  fvx    = vxnew / v0 ;
119  chgx   = MAXI 'ABS' (x - xnew) ;
120  SI (chgx < 0.01) ;
121    info   = CHAI 'It:' &b1 / 5 'Obj:' / 10  psi > 1 'Fvol:' > 4 fvx > 1 'Change:' > 4 chgx > 1 ;
122    MESS info ;
123    QUIT b1 ;
124  FINSI ;
125* préparation de la nouvelle itération
126  x      = xnew ;
127FIN  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. La topologie finale est atteinte après 42 itérations.

../_images/opti_topo_oc.gif

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