************************************************************************ * 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: Method of Moving Asymptotes * * * * Problem: minimize the compliance with: * * a constraint on the global volume fraction * * * * min C(x) = u^T.F * * s.t. G(x) = vf(x)/volfrac - 1 < 0 * * * * | * * | 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. ; ung = MANU 'CHPO' ptg 1 'SCAL' 1. ; ks = kfil * ung ; kfil = NFIL kfil ks ; ** Initialization of the table for the MMA nx = NBEL mesh ; tmma = TABL ; * initial values of the design variables x lx = EXTR x 'VALE' 'SCAL' ; tmma . 'X' = lx ; * bounds for x tmma . 'XMIN' = xmin ; tmma . 'XMAX' = xmax ; * asymptotes tmma . 'LOW' = PROG nx*1. ; tmma . 'UPP' = PROG nx*1. ; * other parameters tmma . 'A0' = 1. ; tmma . 'A' = PROG 0. ; tmma . 'C' = PROG 1.E4 ; tmma . 'D' = PROG 0. ; tmma . 'MOVE' = move ; ** 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)/volfrac - 1 * value g = (vfx / volfrac) - 1. ; * sensitivity (gradient with respect to the physical variable) dg = ve / vtot / volfrac ; * 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 tmma . 'F0VAL' = c ; tmma . 'DF0DX' = EXTR dc 'VALE' 'SCAL' ; tmma . 'FVAL' = PROG g ; tmma . 'DFDX' = TABL ; tmma . 'DFDX' . 1 = EXTR dg 'VALE' 'SCAL' ; MMA tmma ; lxnew = tmma . 'X' ; xnew = MANU 'CHML' mod 'REPA' 'SCAL' lxnew 'TYPE' 'SCALAIRE' 'GRAVITE' ; * 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 ; * summary of the current iteration change = MINI (MAXI 'ABS' (lxnew - (tmma . 'XOLD1'))) (MAXI 'ABS' (lxnew - (tmma . 'XOLD2'))) ; liter = liter ET &b1 ; lobj = lobj ET c ; lvf = lvf ET vfx ; lchange = lchange ET change ; * preparing the next iteration x = xnew ; lx = lxnew ; * ft = 3 --> updatind the slope of the thresholding function SI (ft EGA 3) ; SI ((beta < 512) ET ((loopbeta >EG 50) OU (change 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 ;