Télécharger top_oc.dgibi

Retour à la liste

************************************************************************
* Example of topological optimization                                  *
* Density method + penalization (SIMP)                                 *
*                                                                      *
* This test case is a reproduction of the reference code:              *
* Andreassen, Clausen, Schevenels, Lazarov, & Sigmund (2011)           *
* "Efficient topology optimization in MATLAB using 88 lines of code"   *
* Structural and Multidisciplinary Optimization, 43(1), 1-16           *
* https://doi.org/10.1007/s00158-010-0594-7                            *
*                                                                      *
* Application to a beam under bending in 2D plane stresses             *
*                                                                      *
* Optimization algorithm: Optimality Criteria                          *
*                                                                      *
* Problem: minimize the compliance with a constraint                   *
*          on the volume fraction                                      *
*                                                                      *
*          min   C(x) = u^T.F                                          *
*         s.t.   G(x) = vf(x) = volfrac                                *
*                                                                      *
*       |                                                              *
*       | Force                                                        *
*       v                                                              *
*     O>+------------------------------------------------------+       *
*       | p1                                                   |       *
*       |                                                      |       *
*     O>|                                                      |       *
*       |                                                      |       *
*       |                                                      |       *
*     O>+------------------------------------------------------+ p2    *
*                                                              Ʌ       *
*                                                              O       *
*                                                                      *
************************************************************************
 
 
** 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 ;
 
** Global parameters
itrac    = FAUX ;
OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' 'ECHO' 0 ;
 
** Geometrical parameters (width and height)
l        = 60. ;
h        = 20. ;
e        = 1. ;
 
** Mesh parameters
nelx     = 60 ;
nely     = 20 ;
 
** Material parameters (isotropic elasticity)
e0       = 1. ;
emin     = e0 / 1.E9 ;
nu       = 0.3 ;
 
** Topology optimization parameters
*  penal     : SIMP penalization coefficient
*  volfrac   : minimal volume fraction
*  ft        : = 1 apply filter on the compliance sensitivity field
*              = 2 apply spatial filter on the density field
*              = 3 apply spatial filter on the density field + thresholding by heaviside function
*  rmin      : filter radius (in m)
*  beta      : initial value of the slope of the heaviside function (only for ft = 3)
*              will be doubled every 50 iterations
*  move      : limit of the maximal increment of density
*  changmax  : optimization stop criterion
*  xmin/xmax : min and max density bounds
penal    = 3. ;
volfrac  = 0.5 ;
ft       = 2 ;
rmin     = 0.04 * l ;
beta     = 1 ;
move     = 0.2 ;
changmax = 0.01 ;
xmin     = 0. ;
xmax     = 1. ;
 
** Mesh
p0       = 0. 0. ;
p1       = 0. h ;
ll       = DROI nely p0 p1 ;
mesh     = ll TRAN nelx (l 0.) ;
con      = CONT mesh ;
p2       = con POIN 'PROC' (l 0.) ;
 
** Mechanical model
mod      = MODE mesh 'MECANIQUE' ;
 
** Material properties field with a unit Young modulus
maun     = MATE mod 'YOUN' 1. 'NU' nu 'DIM3' e ;
 
** Boundary conditions
blo      = (BLOQ 'UX' ll) ET (BLOQ 'UY' p2) ;
 
** Load (local force)
f        = FORC (0. -1.) p1 ;
 
** 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. ;
 
** 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 ;
 
** Initial volume fraction
vfx      = (INTG mod xphys maun) / vtot ;
 
** 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. ;
 
** Lets's start a timer
TEMP 'ZERO' ;
 
** Topology optimization loop
liso     = PROG 0. 'PAS' 0.05 1. ;
loopbeta = 0 ;
liter    = PROG ;
lobj     = PROG ;
lvf      = PROG ;
lchange  = PROG ;
REPE b1 500 ;
  loopbeta = loopbeta + 1 ;
* 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 ;
* objective function: compliance = u^T.K.u
* 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 ;
* constraint function: g = vf(x)
* value
  g        = vfx ;
* sensitivity (gradient with respect to the physical variable)
  dg       = ve / vtot ;
* ft = 1  --> filtering the complicance sensitivity
  SI (ft EGA 1) ;
    dc       = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
  FINSI ;
* ft = 2,3 --> updating the sensitivies (to become gradients with respect to the design variable)
  SI (ft EGA 2) ;
    dc       = HFILT dc mod kfil ptg ;
    dg       = HFILT dg mod kfil ptg ;
  FINSI ;
  SI (ft EGA 3) ;
    dx       = (beta * (EXP (-1. * beta * xtilde))) + (EXP (-1. * beta)) ;
    dctilde  = dc * dx ;
    dc       = HFILT dctilde mod kfil ptg ;
    dgtilde  = dg * dx ;
    dg       = HFILT dgtilde mod kfil ptg ;
  FINSI ;
* information about the current topology
  info     = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
  MESS info ;
  SI itrac ;
    def1     = DEFO mesh u 0.05 ;
    TRAC xphys mod con def1 liso 'TITR' info 'NCLK' ;
  FINSI ;
* 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) ;
*   update the physical density (by filtering and thresholding)
    SI (ft EGA 1) ;
      xphys    = xnew ;
    FINSI ;
    SI (ft EGA 2) ;
      xphys    = HFILT xnew mod kfil ptg ;
    FINSI ;
    SI (ft EGA 3) ;
      xtilde   = HFILT xnew mod kfil ptg ;
      xphys    = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
    FINSI ;
*   updating the volume fraction
    vfx      = (INTG mod xphys maun) / vtot ;
    SI (vfx > volfrac) ;
      l1       = lmid ;
    SINON ;
      l2       = lmid ;
    FINSI ;
  FIN b2 ;
* summary of the current iteration
  change   = MAXI 'ABS' (xnew - x) ;
  liter    = liter ET &b1 ;
  lobj     = lobj ET c ;
  lvf      = lvf ET vfx ;
  lchange  = lchange ET change ;
* preparing the next iteration
  x        = xnew ;
* ft = 3 --> updatind the slope of the thresholding function
  SI (ft EGA 3) ;
    SI ((beta < 512) ET ((loopbeta >EG 50) OU (change &lt;EG changmax))) ;
      beta     = 2 * beta ;
      loopbeta = 0 ;
      change   = 1. ;
    FINSI ;
  FINSI ;
* stop criterion
  SI (change < changmax) ;
    info     = CHAI 'It:' &b1 / 5 'Obj:' / 10  c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
    MESS info ;
    QUIT b1 ;
  FINSI ;
FIN  b1 ;
 
** Elapsed time
TEMP 'IMPR' 'SOMM' 'HORL' ;
 
** Plotting the final topology
evobj    = EVOL 'ROUG' 'MANU' 'Iterations' liter 'Compliance' lobj ;
evvf     = EVOL 'ORAN' 'MANU' 'Iterations' liter 'Frac. vol.' lvf ;
evchange = EVOL 'VERT' 'MANU' 'Iterations' liter 'Max. change' lchange ;
SI itrac ;
  info     = CHAI '[Final topology]' ' ' info ;
  TRAC xphys mod con liso 'TITR' info ;
  DESS evobj    'TITR' 'Objective function' ;
  DESS evvf     'TITR' 'Volume fraction' ;
  DESS evchange 'TITR' 'Max. change' ;
FINSI ;
 
** Mesh of the final topology
tab1     = TABL ;
tab1 . 'EPAISSEUR' = e ;
tab1 . 'MODELE'    = mod ;
tab1 . 'TOPOLOGIE' = xphys ;
meshf    = (TOPOSURF tab1) COUL 'GRIS' ;
edge     = ARET meshf ;
SI itrac ;
  TRAC 'FACE' meshf 'ARET' edge 'TITR' 'Mesh of the final topology' ;
FINSI ;
 
FIN ;
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales