Educational article Struct Multidisc Optim   SpringerVerlag  A  line topology optimization code written in Matlab O

Educational article Struct Multidisc Optim SpringerVerlag A line topology optimization code written in Matlab O - Description

Sigmund Abstract The paper presents a compact Matlab im plementation of a topology optimization code for com pliance minimization of statically loaded structures The total number of Matlab input lines is 99 including opti mizer and Finite Element sub ID: 30111 Download Pdf

396K - views

Educational article Struct Multidisc Optim SpringerVerlag A line topology optimization code written in Matlab O

Sigmund Abstract The paper presents a compact Matlab im plementation of a topology optimization code for com pliance minimization of statically loaded structures The total number of Matlab input lines is 99 including opti mizer and Finite Element sub

Similar presentations


Download Pdf

Educational article Struct Multidisc Optim SpringerVerlag A line topology optimization code written in Matlab O




Download Pdf - The PPT/PDF document "Educational article Struct Multidisc Opt..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.



Presentation on theme: "Educational article Struct Multidisc Optim SpringerVerlag A line topology optimization code written in Matlab O"β€” Presentation transcript:


Page 1
Educational article Struct Multidisc Optim 21, 120–127 Springer-Verlag 2001 A 99 line topology optimization code written in Matlab O.Sigmund Abstract The paper presents a compact Matlab im- plementation of a topology optimization code for com- pliance minimization of statically loaded structures. The total number of Matlab input lines is 99 including opti- mizer and Finite Element subroutine. The 99 lines are dividedinto36linesforthemainprogram,12linesforthe Optimality Criteria based optimizer, 16 lines for a mesh- independency filter and 35 lines for the finite

element code. In fact, excluding comment lines and lines associ- ated with output and finite element analysis, it is shown that only )9 Matlab input lines are required for solving a well-posed topology optimization problem. +y adding three additional lines, the program can solve problems with multiple load cases. The code is intended for edu- cational purposes. The complete Matlab code is given in the ,ppendix and can be down-loaded from the web-site http://www.topopt.dtu.dk Keywords topology optimization, education, optimal- ity criteria,world-wideweb,Matlab code Introduction The Matlab

code presented in this paper is intended for engineering education. -tudents and newcomers to the field of topology optimization can down-load the code from the web-page http://www.topopt.dtu.dk The code may be used in courses in structural optimiza- tion where students may be assigned to do extensions such as multiple load-cases, alternative mesh-independ- ency schemes,passiveareas,etc. ,notherpossibility is to usetheprogramtodevelopstudents.intuitionforoptimal design. ,dvanced students may be askedto guess the op- timal topology for given boundary condition and volume Received October

22, 1999 O.Sigmund Department of Solid Mechanics, Building 404, Technical Uni- versityofDenmark,DK-2800(yngby,Denmark e-mail: sigmund@fam.dtu.dk fraction and then the programshows the correctoptimal topologyforcomparison. Intheliterature,onecanfindamultitudeofapproaches for the solving of topologyoptimization problems. In the original paper +ends0e and 1ikuchi 219334 a so-called microstructure or homogenization based approach was used, based on studies of existence of solutions. Thehomogenizationbasedapproachhasbeenadopted in many papers but has the disadvantage that the deter- mination

and evaluation of optimal microstructures and their orientations is cumbersome if not unresolved 2for noncompliance problems4 and furthermore, the resulting structures cannot be built since no definite length-scale is associated with the microstructures. 5owever, the ho- mogenization approach to topology optimization is still important in the sense that it can provide bounds on the theoreticalperformanceofstructures. ,n alternative approach to topology optimization is the so-called 6power-law approach7 or -IM8 approach 2-olid Isotropic Material with 8enalization4 2+ends0e 19399:hou and

;ozvany 19919Mlejnek 19924. 5ere, ma- terial properties are assumed constant within each elem- entusedtodiscretizethedesigndomainandthevariables are the element relative densities. The material proper- ties are modelled as the relative material density raised to some power times the material properties of solid ma- terial. This approach has been criticized since it was ar- gued that no physical material exists with properties de- scribedbythepower-lawinterpolation.5owever,arecent paper by +ends0e and -igmund 219994 proved that the power-law approach is physically permissible as long as simple

conditions on the power are satisfied 2e.g. for 8oisson.s ratio equal to 4. To ensure existence of so- lutions, the power-law approach must be combined with a perimeter constraint, a gradient constraint or with fil- tering techniques 2see -igmund and 8etersson 1993, for an overview4. The power-law approach to topology op- timization has been applied to problems with multiple constraints,multiple physicsand multiple materials. =hereas the solution of the above mentioned ap- proaches is based on mathematical programming tech- niques and continuous design variables, a number of pa-

pers have appeared on solving the topology optimization problem as an integer problem. +eckers 219994 success-
Page 2
121 fully solved large-scale compliance minimization prob- lems using a dual-approach but other approaches based on genetic algorithms or other semi-random approaches require thousands of function evaluations even for small number ofelements andmustbe consideredimpractical. ,part from above mentioned approaches, which all solve well defined problems 2e.g. minimization of com- pliance4 a number of heuristic or intuition based ap- proaches have been shown to

decrease compliance or other objective functions. ,mong these methods are so- called evolutionary design methods 2see e.g. Xie and -teven 19979 +aumgartner etal. 19924. ,part from be- ing very easy to understand and implement 2at least for the compliance minimization case4, the main moti- vation for the evolutionary approaches seems to be that mathematically based or continuous variable approaches 6involve some complex calculus operations and mathe- matical programming7 2citation from Li etal. 19994 and they contain 6mathematical methods of some complex- ity7 2citation from :hao etal. 19934

whereas the evo- lutionary approach 6takes advantage of powerful com- puting technology and intuitive concepts of evolution processes in nature7 2citation from Li etal. 19994. Two things can be arguedagainstthis. First,the evolutionary approaches become complicated themselves, once more complex objectives than compliance minimization are consideredandsecond,asshowninthispaper,the6math- ematically based7 approaches for compliance minimiza- tion are simple to implement as well and are compu- tationally equally eAcient. Furthermore, mathematical programming based methods can easily be extended to

other non-compliance objectives such as non-self-adjoint andmultiphysicsproblemsandtoproblemswithmultiple constraints2e.g. -igmund19994.Extensions ofthe evolu- tionaryapproachtosuchcasesseemmorequestionable. The complete Matlab code is given in the ,ppendix. The remainder of the paper consists of definition and discussion of the optimization problem 2-ect. 24, com- ments about the Matlab implementation 2-ect. 34 fol- lowedbyadiscussionofextensions2-ect.)4andaconclu- sion2-ect. 54. The topology optimization problem ,numberofsimplificationsareintroducedtosimplifythe Matlab code.

First, the design domain is assumed to be rectangular and discretized by square finite elements. In this way, the numbering of elements and nodes is simple 2column by column starting in the upper left corner4 and the aspect ratio of the structure is given by the ratio of elements in the horizontal 2 nelx 4 and the vertical direc- tion 2 nely 4. Names in type-writer style refer to Matlab variable names that di-er from the obvious .see the Matlab code in the /ppendix1 ,topologyoptimizationproblembasedonthepower- law approach, where the objective is to minimize compli- ancecanbe writtenas

min 4C KU =1 subjecttoB KU B0 min 214 where and are the global displacement and force vectors, respectively, is the global stiEness matrix, and are the element displacement vector and stiEness matrix, respectively, is the vector of design variables, min is a vector of minimum relative densities 2non-zero to avoid singularity4, 2C nelx nely 4isthenumber of elements used to discretize the design domain, is the penalization power 2typically C34, 4and is the material volume and design domain volume, respectively and volfrac 4is the prescribedvolumefraction. The optimization problem 214 could be

solved using several diEerent approaches such as Optimality Criteria 2OC4 methods, -equential Linear 8rogramming 2-L84 methodsortheMethodofMoving,symptotes2MM,by -vanberg 19374 and others. For simplicity, we will here useastandardOC-method. Following+ends0e219954aheuristicupdatingscheme forthedesignvariablescanbeformulatedas new max2 min ,x if max2 min ,x if max2 min ,x min21 ,x min21 ,x if min21 ,x 224 where move 4 is a positive move-limit, 2C 1 24 is anumericaldampingcoeAcientand is found fromthe optimalityconditionas ∂c ∂x ∂V ∂x 234 where is a Lagrangian multiplier

that can be found by abi-sectioningalgorithm. The sensitivityofthe objectivefunction isfound as ∂c ∂x 2)4
Page 3
122 Formoredetails onthederivationandimplementationof the optimality criteria method, the reader is referred to theliterature2e.g.+ends0e19954. In order to ensure existence of solutions to the top- ology optimization problem 214, some sort of restriction on the resulting design must be introduced 2see -igmund and 8etersson 1993, for an overview4. 5ere we use a fil- tering technique 2-igmund 199), 19974. It must be em- phasized that this filter has

not yet been proven to en- sure existence of solutions, but numerous applications by the author have proven the the filter produces mesh- independent designsin practice. Themesh-independencyfilterworksbymodifyingthe elementsensitivities asfollowsB ∂c ∂x =1 =1 ∂c ∂x 254 Theconvolutionoperator2weightfactor4 iswrittenas min dist2 e,f dist2 e,f min ,e C1 ,...,N, 264 wheretheoperatordist2 e,f 4isdefinedasthedistancebe- tween centre of element and centre of element .The convolution operator is zero outside the filter area. The convolution operator

decays linearly with the dis- tance from element . Instead of the original sensitivities 2)4,themodifiedsensitivities254areusedintheOptimal- ityCriteriaupdate 234. Matlab implementation TheMatlabcode2seethe,ppendix4,isbuiltupasastan- dard topology optimization code. The main program is called fromthe Matlabprompt bythe line top(nelx,nely,volfrac,penal,rmin) where nelx and nely are the number of elements in the horizontal and vertical directions, respectively, volfrac is the volume fraction, penal is the penalization power and rmin isthefiltersize2dividedbyelementsize4.Other

variables as well as boundary conditions are defined in the Matlab code itself and can be edited if needed. For eachiterationinthetopologyoptimizationloop,thecode generates a picture of the current density distribution. Figure1showstheresultingdensitydistributionobtained by the code given in the ,ppendix called with the input line top(60,20,0.5,3.0,1.5)         Fig.1 Topology optimization of the MBB-beam. Top: full design domain, middle: half design domain with symmetry boundary conditions and bottom: resulting topology opti- mized beam .both halves1

Thedefaultboundaryconditionscorrespondtohalfofthe 6M++-beam7 2Fig. 14. The load is applied vertically in the upper left corner and there is symmetric boundary conditions along the left edge and the structure is sup- portedhorizontallyinthelowerrightcorner. Important details of the Matlab code are discussed in the followingsubsections. 3.1 Main program (lines 1–36) The main program2lines 1–364starts by distributing the material evenly in the design domain 2line )4. ,fter some other initializations, the main loop starts with a call to theFiniteElementsubroutine2line124whichreturnsthe

displacement vector . -ince the element stiEness matrix for solid material is the same for all elements, the elem- ent stiEness matrix subroutine is called only once 2line 1)4. Following this, a loop over all elements 2lines 16 2)4 determines objective function and sensitivities 2)4. The variables n1 and n2 denote upper left and right element node numbers in global node numbers and are used to extract the element displacement vector Ue from the global displacement vector . The sensitivity analy- sis is followed by a call to the mesh-independency filter 2line 264 and the Optimality

Criteria optimizer 2line 234. The current compliance as well as other parameters are printed by lines 30–33 and the resulting density distri- bution is plotted 2line 354. The main loop is terminated if the change in design variables 2 change determined in line304is lessthan1percent . Otherwiseabovestepsare repeated. this is a rather 3sloppy4 convergence criterion and could be decreased if needed
Page 4
123 3.2 Optimality criteria based optimizer (lines 37–48) The updated design variables are found by the opti- mizer 2lines 37–)34. 1nowing that the material volume sum(sum(xnew)) 4 is

a monotonously decreasing func- tion of the Lagrange multiplier 2 lag 4, the value of the Lagrangianmultiplierthatsatisfiesthevolumeconstraint can be found by a bi-sectioning algorithm 2lines )0-)34. The bi-sectioning algorithm is initialized by guessing alower l1 andanupper l2 bound for the Lagrangian multiplier 2line 394. The interval which bounds the La- grangian multiplier is repeatedly halved until its size is lessthan the convergencecriteria2line )04. 3.3 Mesh-independency (ltering (lines 49–64) Lines )9–6) represent the Matlab implementation of 254. Iote that not all elements in

the design domain are searched in order to find the elements that lie within the radius rmin but only those within a square with side lengths two times round(rmin) around the considered element. +y selecting rmin less than oneinthe callofthe routine, the filtered sensitivities will be equal to the ori- ginalsensitivities makingthe filterinactive. 3.4 )inite element code (lines 65–99) The finite element code is written in lines 65–99. Iote that the solver makes use of the sparse option in Mat- lab. The global stiEness matrix is formed by a loop over all elements 2lines

70–774. ,s was the case in the main program,variables n1 and n2 denote upper left and right element node numbers in global node numbers and are used to insert the element stiEness matrix at the right placesinthe globalstiEness matrix. ,s mentioned before, both nodes and elements are numbered column wise from left to right. Furthermore, eachnodehastwodegreesoffreedom2horizontalandver- tical4, thus the command F(2,1)=-1. 2line 794 applies averticalunitforceforceinthe upper leftcorner. -upports are implemented by eliminating fixed de- greesoffreedomfromthelinearequations.Matlabcando this

veryelegantlywith the line 8( U(freedofs,:) = )(freedofs,freedofs) * F(freedofs,:)+ where freedofs indicate the degrees of freedom which are unconstrained. Mostly, it is easier to define the de- greesoffreedomthatarefixed2 fixeddofs 4thereafterthe freedofs arefoundautomaticallyusingtheMatlaboper- ator setdiff which finds the free degrees of freedoms as         Fig.2 Topology optimization of a cantilever beam. (eft: de- sign domain and right: topology optimized beam thediEerencebetweenalldegreesoffreedomandthefixed degreesoffreedom2line324. The element

stiEness matrix is calculated in lines 36 99. The 3 by 3 matrix for a square bi-linear )-node elem- ent was determined analytically using a symbolic manip- ulation software. The Young.s modulus and the 8ois- son.sratio nu canbe alteredinlines33and39. Extensions The Matlab code given in the ,ppendix solves the prob- lem of optimizing the material distribution in the M++- beam 2Fig. 14 such that its compliance is minimized. , number of extensions and changes in the algorithm can be thought of, a few of which are mentioned in the following. 4.1 Other boundary conditions It is very simple to

change boundary conditions and sup- portconditionsinordertosolveotheroptimizationprob- lems.Inordertosolvetheshortcantileverexampleshown inFig.2,onlylines79and30mustbechangedto -. F(2/(nelx01)/(nely01),1) = -1+ 80 fixeddofs = 11:2/(nely01)2+ =ith these changes, the input line for the case shown in Fig.2is top(32,20,0.(,3.0,1.2) 4.2 Multiple load cases It is also very simple to extend the algorithm to account formultipleloadcases.Infact,thiscanbedonebyadding only three additional lines and making minor changes to another) lines.
Page 5
12) In the case of two load cases, force and

displacement vectors must be defined as two-column vectors which meansthat line 69is changedto 6. F = sparse(2/(nely01)/(nelx01),2)+ U = sparse(2/(nely01)/(nelx01),2)+ The objective function is now the sum of two compli- ances,i.e. 4C =1 KU 274 thus lines 20–22aresubstitutedwiththelines 1.3 dc(ely,elx) = 0.+ 1.c for i = 1:2 20 Ue = U(12/n1-1+2/n1+ 2/n2-1+2/n2+ 2/n201+2/n202+2/n101+2/n1022,i)+ 21c = c 0 x(ely,elx)4penal/Ue5/),/Ue+ 22 dc(ely,elx) = dc(ely,elx) - penal/x(ely,elx)4(penal-1)/Ue5/),/Ue+ 223 end To solve the two-load problem indicated in Fig. 3, a unit upward load in the

top-right corner is added to line 79, whichthenbecomes -. F(2/(nelx01)/(nely01),1) = -1.+ F(2/(nelx)/(nely01)02,2) = 1.+ The input line forFig. 3is top(30,30,0.(,3.0,1.2). 4.3 Passi-e elements In some cases, some of the elements may be required to takethe minimum densityvalue 2e.g.a holeforapipe4. ,n nely nelx array passive with zerosatelements free to change and ones at elements fixed to be zero can bedefinedinthemainprogramandtransferredtotheOC subroutine2adding passive tothecallinlines23and334. The addedline (23 xnew(find(passive)) = 0.001+ in the OC subroutine looks for passive

elements and sets their densityequal tothe minimum density20.0014. 22         Fig.3 Topology optimization of a cantilever beam with two load-cases. (eft: design domain, middle: topology optimized beam using oneload case and right: topology optimized beam using two load cases Fig.4 Topology optimization of a cantilever beam with a 5xed hole. (eft: design domain and right: topology opti- mized beam Figure ) shows the resulting structure obtained with the input top((5,30,0.5,3.0,1.5), when the following 10 lines were added to the main program 2after line )4 in order to find

passive elem- ents within a circle with radius nely/3. and center nely/2. nelx/3. for ely = 1:nely for elx = 1:nelx if s6rt((ely-nely/2.)420(elx-nelx/3.)42) 7 nely/3. passive(ely,elx) = 1+ x(ely,elx) = 0.001+ else passive(ely,elx) = 0+ end end end 4.4 Alternati-e optimizer ,dmittedly, the optimality criteria based optimizer im- plemented here is only good for a single constraintand it is based on a heuristic fixed point type updating scheme. Inordertoinstallabetteroptimizer,onecanobtain2free
Page 6
125 of charge for academic purposes4 the Matlab version of the MM,-algorithm

2-vanberg 19374 from 1rister -van- berg, 1T5, -weden. The MM, code is called with the followinginput line mmasu3(89:U;-varia3les, ... , where the total number of inputKoutput variables is 20, including objective function, constraints, old and new densities, etc. Implementing the MM,-optimizeris fairly simple, but requires the definition of several auxiliary variables. 5owever,it allows for the solvingof more com- plexdesignproblemswithmorethanoneconstraint.The Matlab optimizer will solve the standard topology op- timization problem using less iterations at the cost of aslightly

increasedC8L-time periteration. 4.5 Other extensions Extensions to three dimensions should be straight for- ward whereas more complex problems such as compliant mechanism design 2-igmund 19974 requires the imple- mentationoftheMM,optimizerandthedefinitionofex- tra constraints. The simplicity of the Matlab commands allowforeasyextensionsofthe graphicaloutput, interac- tive inputetc. .onclusions This paper has presented a very simple implementation of a mathematicalprogrammingbase topologyoptimiza- tion algorithm. The code is implemented using only 99 Matlab input lines and includes

optimizer, mesh-indepe- ndency filteringand Finite Elementcode. The Matlab code can be down-loaded from the web- page http://www.topopt.dtu.dk andisintendedfored- ucational purposes. The code can easily be extended to include multi load problems and the definition of passive areas. ;unning the code in Matlab is rather slow compared to a Fortran implementation of the same code which can be tested at the web-site http://www.topopt.dtu.dk 5owever,an add-onpackageto Matlab 2M,TL,+ Com- piler4 allows for the generation of more eAcient C-code that can be optimized for run-time 2this

option, how- ever, has not been tested by the author4. It should be noted that speed can be gained by modifying the Mat- lab code itself, howeverthe speed is gained on the cost of simplicity of the program. The modification is suggested by ,ndreas ;ietz from LinkM oping Lniversity who uses sparsityoptionsintheassemblyoftheglobalstiEnessma- trix.Thereadermaydown-loadhiscodeattheweb-pageB http://www.mekanik.ikp.liu.se/andridiv/matla3/ theory.html The code was intentionally kept compact in order to keep the total number of lines below 100. If users of the code should find ways to

further compactify or simplify the code, the author would be happy to re- ceive suggested modifications that can be implemented in the public domain code 2the author.s e-mail address is sigmund=fam.dtu.dk 4. -ince its first publication on the =orld =ide =eb in October 1999, the Matlab code has been down-loaded more than 500 times by diEerent users 2as of ,ugust 20004. ,mong other positive feedbacks, several profes- sors reported that they have used the code in courses on structuraloptimizationandhavelettheirstudentsimple- ment alternative boundary conditions and multiple load

cases. Acknowledgements This work was supported by the Dan- ish Technical 6esearch 7ouncil through the T8O69Talent- programme: Design of Micro:lectroMechanical Systems .M:MS1. The author would also like to thank Thomas Buhl, Technical University of Denmark, for his inputs to an earlier version of the code. References Baumgartner, /., 8arzheim, (.; Mattheck, 7. 1992: SKO .Soft Kill Option1: The biological way to 5nd an optimum structure topology. Int. J. Fatigue 14 , 387–393 Beckers,M.1999: Topology optimization usingadualmethod with discrete variables. Struct. Optim. 17 , 14–24 Bendse,

M.P. 1989: Optimal shape design as a material dis- tribution problem. Struct. Optim. , 193–202 Bendse, M.P. 1995: Optimization of structural topology, shape and material . Berlin, 8eidelberg, New Aork: Springer Bendse, M.P.; Kikuchi, N. 1988: Benerating optimal topolo- giesinoptimaldesign usingahomogenization method. Comp. Meth. Appl. Mech. Engrg. 71 , 197–224 Bendse, M.P.; Sigmund, O. 1999: Material interpolations in topology optimization. Arch. Appl. Mech. 69 , C35–C54 (i, D.;Steven,B.P.; Eie, A.M. 1999: On equivalencebetween stress criterion and sti-ness criterion in

evolutionary struc- tural optimization. Struct. Optim. 18 , C7–73 MleGnek, 8.P. 1992: Some aspects of the genesis of structures. Struct. Optim. , C4–C9 Sigmund, O. 1994: Design of material structures using top- ology optimization . Ph.D. Thesis, Department of Solid Me- chanics, Technical University of Denmark Sigmund, O. 1997: On the design of compliant mechanisms using topology optimization. Mech. Struct. Mach. 25 , 495 52C
Page 7
126 Sigmund, O. 1999: Topology optimization of multi-physics, multi-material structures. Proc. WCSMO-3 .heldinBu-alo, NA1 Sigmund, O.; Petersson, H. 1998:

Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct. Optim. 16 , C8–75 Svanberg,K.1987:Themethodofmovingasymptotes–anew method for structural optimization. Int. J. ,umer. Meth. En- grg. 24 , 359–373 Eie, A.M.; Steven, B.P. 1997: Evolutionary structural opti- mization . Berlin, 8eidelberg, New Aork: Springer Ihao, 7.; 8ornby, P.; Steven, B.P.; Eie, A.M. 1998: / gener- alized evolutionary method for numerical topology optimiza- tion of structures under static loading conditions. Struct. Op- tim. 15 ,

251–2C0 Ihou, M.; 6ozvany, B.J.N. 1991: The 7O7 algorithm, part JJ:Topological, geometryandgeneralized shapeoptimization. Comp. Meth. Appl. Mech. Engrng. 89 , 197–224 Appendix –Matlab code 1 %%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, OCTOBER 1999 %%% 2 function top(nelx,nely,volfrac,penal,rmin/0 1 % INITIALIZE 2 x(13nely,13nelx/ 4 volfrac0 5 loop 4 60 7 change 4 1.0 9 % START ITERATION : while change > 6.61 9 loop 4 loop = 10 16 xold 4 x0 11 % >E-ANALYSIS 12 @UA4>E(nelx,nely,x,penal/0 11 % OBBECTICE >UNCTION AND SENSITICITY ANALYSIS 12 @DEA 4 lk0 15 c 4 6.0 17 for ely 4 13nely

19 for elx 4 13nelx 1: n1 4 (nely=1/E(elx-1/=ely0 19 n2 4 (nely=1/E elx =ely0 26 Ue 4 U(@2En1-102En10 2En2-102En20 2En2=10 2En2=20 2En1=102En1=2A,1/0 21 c 4 c = x(ely,elx/FpenalEUeGEDEEUe0 22 dc(ely,elx/ 4 -penalEx(ely,elx/F(penal-1/E UeGEDEEUe0 21 end 22 end 25 % >ILTERING O> SENSITICITIES 27 @dcA 4 check(nelx,nely,rmin,x,dc/0 29 % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD 2: @xA 4 OC(nelx,nely,x,volfrac,dc/0 29 % PRINT RESULTS 16 change 4 max(max(aIs(x-xold///0 11 disp(@G It.3 G sprintf(G%2iG,loop/ G OIJ.3 G sprintf(G%16.2fG,c/ ... 12 G Col.3 G sprintf(G%7.1fG,sum(sum(x//K (nelxEnely//

... 11 G ch.3 G sprintf(G%7.1fG,change /A/ 12 % PLOT DENSITIES 15 colormap(gray/0 imagesc(-x/0 axis eLual0 axis tight0 axis off0pause(1e-7/0 17 end 19 %%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%% 1: function @xnewA4OC(nelx,nely,x,volfrac,dc/ 19 l1 4 60 l2 4 1666660 move 4 6.20 26 while (l2-l1 > 1e-2/ 21 lmid 4 6.5E(l2=l1/0 22 xnew 4 max(6.661,max(x-move,min(1.,min(x=move,x. EsLrt(-dc.Klmid/////0 21 if sum(sum(xnew// - volfracEnelxEnely > 60 22 l1 4 lmid0 25 else 27 l2 4 lmid0 29 end 2: end 29 %%%%%%%%%% MESH-INDEPENDENCY >ILTER %%%%%%%%%%% 56 function @dcnA4check(nelx,nely,rmin,x,dc/ 51

dcn4zeros(nely,nelx/0 52 for i 4 13nelx 51 for J 4 13nely 52 sum46.60 55 for k 4 max(i-round(rmin/,1/3 min(i=round(rmin/,nelx/ 57 for l 4 max(J-round(rmin/,1/3 min(J=round(rmin/, nely/ 59 fac 4 rmin-sLrt((i-k/F2=(J-l/F2/0 5: sum 4 sum=max(6,fac/0 59 dcn(J,i/ 4 dcn(J,i/ = max(6,fac/Ex(l,k/ Edc(l,k/0 76 end 71 end 72 dcn(J,i/ 4 dcn(J,i/K(x(J,i/Esum/0 71 end 72 end 75 %%%%%%%%%% >E-ANALYSIS %%%%%%%%%%%% 77 function @UA4>E(nelx,nely,x,penal/ 79 @DEA 4 lk0 7: D 4 sparse(2E(nelx=1/E(nely=1/, 2E(nelx=1/E (nely=1//0 79 > 4 sparse(2E(nely=1/E(nelx=1/,1/0 U 4 sparse(2E(nely=1/E(nelx=1/,1/0 96 for ely 4

13nely 91 for elx 4 13nelx 92 n1 4 (nely=1/E(elx-1/=ely0 91 n2 4 (nely=1/E elx =ely0 92 edof 4 @2En1-10 2En10 2En2-10 2En20 2En2=10 2En2=202En1=10 2En1=2A0 95 D(edof,edof/ 4 D(edof,edof/ = x(ely,elx/FpenalEDE0 97 end 99 end 9: % DE>INE LOADSAND SUPPORTS(HAL> MBB-BEAM/ 99 >(2,1/ 4 -10 :6 fixeddofs 4 union(@13232E(nely=1/A, @2E(nelx=1/E(nely=1/A/0 :1 alldofs 4 @132E(nely=1/E(nelx=1/A0 :2 freedofs 4 setdiff(alldofs,fixeddofs/0 :1 % SOLCING
Page 8
127 :2 U(freedofs,3/ 4 D(freedofs,freedofs/ N >(freedofs,3/0 :5 U(fixeddofs,3/4 60 :7 %%%%%%%%%% ELEMENT STI>>NESS MATRIO %%%%%%% :9 function

@DEA4lk :: E 4 1.0 :9 nu 4 6.10 96 k4@ 1K2-nuK7 1K:=nuK: -1K2-nuK12 -1K:=1EnuK: ... 91 -1K2=nuK12 -1K:-nuK: nuK7 1K:-1EnuK:A0 92 DE 4 EK(1-nuF2/E @ k(1/ k(2/ k(1/ k(2/ k(5/ k(7/ k(9/ k(:/ 91 k(2/ k(1/ k(:/ k(9/ k(7/ k(5/ k(2/ k(1/ 92 k(1/ k(:/ k(1/ k(7/ k(9/ k(2/ k(5/ k(2/ 95 k(2/ k(9/ k(7/ k(1/ k(:/ k(1/ k(2/ k(5/ 97 k(5/ k(7/ k(9/ k(:/ k(1/ k(2/ k(1/ k(2/ 99 k(7/ k(5/ k(2/ k(1/ k(2/ k(1/ k(:/ k(9/ 9: k(9/ k(2/ k(5/ k(2/ k(1/ k(:/ k(1/ k(7/ 99 k(:/ k(1/ k(2/ k(5/ k(2/ k(9/ k(7/ k(1/A0