Aconcurrentﬁbreorientationandtopologyoptimisationframeworkfor3D-printedcarbonﬁbre- reinforced composites José Humberto S. Almeida Jr.,Bruno G. Christoﬀ,Volnei Tita,Luc St-Pierre Isotropic OrthotropicStructural stiffness+250% 1. Design 2. Optimisation 3. Manufacturing 4. Testing Aconcurrentﬁbreorientationandtopologyoptimisationframeworkfor 3D-printed carbon ﬁbre-reinforced composites José Humberto S. Almeida Jr.a,b,<, Bruno G. Christoﬀc, Volnei Titac,dand Luc St-Pierreb aSchool of Mechanical and Aerospace Engineering, Queen’s University Belfast, Belfast, UK bDepartment of Mechanical Engineering, Aalto University, Espoo, Finland cDepartment of Aeronautical Engineering, São Carlos School of Engineering, University of São Paulo, São Carlos, SP, Brazil dDepartment of Mechanical Engineering, Faculty of Engineering of the University of Porto, Porto, Portugal ARTICLE INFO Keywords : Topology optimisation Material optimisation Concurrent optimisation Additive manufacturingABSTRACT This work proposes a novel framework able to optimise both topology and ﬁbre angle concomitantly to maximise the stiﬀness of a structure. Two diﬀerent materials are considered, one with isotropic properties (nylon) and another one with orthotropic properties (onyx, which is nylon reinforced with choppedcarbonﬁbres). Theframeworkoptimises,inthesameparticularsub–step,ﬁrstthetopology, and second, the ﬁbre angle at every element throughout the design domain. For the isotropic mate- rial,onlytopologyoptimisationtakesplace,whereasfortheorthotropicsolid,bothtopologyandﬁbre orientation are considered. The objective function is to minimise compliance and three admissible volumes: 30%,40%,and50%. Threeclassicalbenchmarkcasesareconsidered: acantileverbeam,as well as 3-point and 4-point bending. The optimum topologies are further treated and manufactured using the fused ﬁlament fabrication (FFF) 3D printing method. Key results reveal that the absolute stiﬀness,density–normalisedandvolume–normalisedstiﬀnessvalueswithineachadmissiblevolume are higher for onyx than for nylon, which proves the eﬃciency of the proposed concurrent optimisa- tion framework. Moreover, although the objective function was to minimise compliance, it was also eﬀectivetoimprovethestrengthofallparts. Theexcellentqualityandgeometrictoleranceofthe3D printed parts are also worth mentioning. Highlights •A topology optimisation framework for isotropic and orthotropic materials. •Consideration of ﬁbre orientation into the optimisation framework. •Additive manufacturing of chopped carbon ﬁbre reinforced nylon composites. •Experimental validation: manufacturing and testing. 1. Introduction One of the most fundamental engineering challenge is how to design a structure to be as light as possible with- out sacriﬁcing its mechanical performances. According to SigmundandMaute[1],thiscanbeachievedwithtopology optimisation (TO). Its basic idea relies on repeated analysis and design update steps usually guided by a gradient-based computation. The ﬁrst attempt on TO was carried out by BendsøeandKikuchi[2]withtheaimtoproposeanalterna- tivemethodtoshapeoptimisationapproachabletoyieldboth the optimum topology and the optimal shape of a structure. Later on, Bendsøe [3] introduced several ways of removing the discrete nature of the problem by introducing a density functionasacontinuousdesignvariablewithintheoptimisa- tionproblem. Infact,BendsøeandKikuchi[2],Bendsøe[3] andRozvanyetal. [4]proposedthemostdisseminatedmath- ematical approach for TO, the well-known Solid Isotropic <Corresponding Author humberto.almeida@qub.ac.uk (J.H.S. Almeida Jr.); brunogch@gmail.com (B.G. Christoﬀ); voltita@sc.usp.br (V. Tita); luc.st-pierre@aalto.fi (L. St-Pierre) ORCID(s): 0000-0002-9408-7674 (J.H.S. Almeida Jr.); 0000-0003-2148-4596 (B.G. Christoﬀ); 0000-0002-8199-1162 (V. Tita); 0000-0003-3857-5398 (L. St-Pierre)MaterialwithPenalisation(SIMP)method[5]. Thismethod ﬁndstheoptimalmaterialdistributionwithinaparticularde- sign domain, load cases, boundary conditions, manufactur- ingconstraints,andperformancerequirements[6,7]. Other less explored approaches include, for instance, the level set [8], evolutionary structural optimisation [9], and moving morphable component methods [10]. Inthelast30years,severalTOapproacheshavebeenex- plored, mostly for isotropic materials [11]. This is mainly due to two reasons: i) the high maturity level of ap- proachesforisotropicmaterials(suchasSIMP)andii)well- established manufacturing techniques for complex metallic shapes. However, the continuous demand for lightweight structures has led to the increase use of anisotropic carbon ﬁbre-reinforced polymer (CFRP) composites, especially in high-performance aerospace, aeronautical, and automotive components [12]. From the manufacturing point of view, complex CFRP shapes can now be produced by additive manufacturing (AM) techniques [13, 14]. Among them, fused deposition modelling (FDM), or fused ﬁlament fab- rication (FFF), is one of the most disseminated AM tech- niques, in which a polymer ﬁlament, often reinforced with ﬁbres,isextrudedanddepositedonabuildplate[15,16,17]. JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 1 of 15 Concurrent topology and ﬁbre angle optimisation of composites According to the systematic review of Sanei and Popescu [18], most studies on FFF for CFRP parts were published in the last ﬁve years only [19]. For instance, Chen and Ye [20, 21] 3D printed composite parts with carbon ﬁbre and nylon materials using topology optimisation relying on the classical SIMP method [20] and to 3D print structures with negative Poisson’s ratio [21]. Sugiyama et al. [22] 3D printed optimised composites with carbon ﬁbres placed along the principal stress directions. However, there is still an open question: how to arrange the reinforcing ﬁbres and optimise the part topology? Aim- ing at taking full advantage of both material anisotropy and exploiting the powerful manufacturing capability of FFF method,anadequateanswerisbycombiningTOtooptimise the material distribution with structural optimisation to ad- justthelocalﬁbreorientation. Ontheonehand,thematerial distributioncanbeoptimisedusingSIMP;ontheotherhand, theﬁbreanglecanbelocallyoptimisedwithgradient-based algorithms, in which the ﬁbre angle can be treated as a de- sign variable [23]. Thereareafewstudiesdealingwithbothtopologyandﬁ- breangleoptimisationssimultaneously. Forexample,Leeet al. [24] developed a TO method for optimising both mate- riallayoutandﬁberorientationinfunctionallygradedﬁber- reinforced composites, in which the ﬁbre angle are consid- eredasdiscretevariables. Tongetal. [25]builtasequential optimisationmethodconsideringbothﬁbreangleandtopol- ogy for constant-stiﬀness laminated plates using lamination parametersasdesignvariables. Thestiﬀnessforashortcan- tilever beam and the ﬂexibility for a compliant inverter in- creased by 6.5% and 4.2%, respectively. Jiang et al. [26] proposed a TO approach for continuous ﬁber angle optimi- sation, which computes the best layout and orientation of ﬁber reinforcement for AM structures. They report a mini- mum compliance 63% lower than the baseline by selecting a diﬀerent print orientation, in which the ﬁber orientation follows the outer contour of the dense material region for each layer. Nomura et al. [27] developed a TO framework for designing both topology and orientation distributions of composite materials simultaneously. However, the optimi- sation results have some areas which violate the conditions for tensor invariants. Papapetrou et al. [28] developed an optimisationframeworkforbothtopologyandﬁbrepathsto create variable-stiﬀness designs. In general, the optimised part is stiﬀer than the baseline. Yan et al. [29] proposed a concurrenthierarchicaloptimisationmethodologyconsider- ingthesimultaneousoptimisationofstructuraltopologyand orthotropic material orientation. They showed by means of numerical examples that optimising both topology and ﬁ- bre angle might decrease compliance. Among these stud- ies, only Jiang et al. [26] manufactured optimum beams us- ing carbon-ﬁber-reinforced polylactic acid (PLA) compos- ites. However, only one type of specimen and one volume was3D-printed. Nevertheless,alltheothermentionedstud- ies are merely computational. Chen and Ye [20] developed a procedure combining topological design and ﬁbre place- ment paths based on average load transmission trajectoriesfor3D-printedCFRPpartsusingtheclassicalSIMPmethod for optimising the topology. After the non–extensive state–of-the–art on research un- dertaken on the topic, the following gaps have been identi- ﬁed, which underlie the realisation of this study: •no studies are considering the anisotropy of carbon ﬁ- bres in the topology optimisation; •there is need for a framework that can simultaneously optimise both topology and ﬁbre angle; •no investigations are considering chopped carbon ﬁbre reinforced thermoplastic composites; and •there is a lack of comprehensive investigations consid- ering the whole design process: from optimisation to manufacturing and testing. In this context, this work proposes a robust concurrent topology and ﬁbre angle optimisation framework for 3D- printed composites using dedicated algorithms considering the orthotropy of the composite materials in contrast to the classicalSIMPapproach. Thetopologyproblemissolvedas a constrained problem whereas the ﬁbre angle is treated as anunconstrainedoptimisationproblem. Threeﬁnalvolumes are selected (30%, 40%, and 50%), and the objective func- tioniscomplianceminimisation. Threebenchmarkcasesare consideredtoevaluatethenovelapproachhereindeveloped: i)cantileverbeam,ii)three-pointbending,andiii)four-point bending(orMMBbeam). Twomaterialspropertiesarecon- sidered: isotropic and orthotropic. The isotropic material consists of nylon (Polyamide 6 – nylon), whereas the trans- verselyisotropicmaterialisachoppedcarbon-ﬁbreﬁlledny- lon, called onyx. After the numerical optimisation is ﬁn- ished, the optimal 3-point and 4-point bending beams are 3D printed, with both nylon and onyx materials, and tested under the same boundary conditions (BC) used in the opti- misation. 2. The formulation Let be a domain in R2with boundary ) (Figure 1). Consider that the Dirichlet and Neumann boundary condi- tions are applied on ) uand) f, respectively, where ) =) uiä) fi;) uiã) fi=Éi=1;2;(1) such that ui=ug ion) uii=1;2; (2) whereuiisthei-thcomponentofthedisplacementand ug iis known, and ijnj=tion) fii=1;2; (3) JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 2 of 15 Concurrent topology and ﬁbre angle optimisation of composites Ω𝑓 𝜕Ω void 𝜌=0 𝑥𝑦 solid 𝜌=1𝜃1 𝜃3𝜃2 Figure 1: Schematic representation of both topology and ﬁbre orientation within a given domain. whereis the stress tensor, nis a normal unit vector and ti is the traction vector. Wewanttoﬁndtheoptimalmaterialdistributionthatmin- imises compliance when subjected to boundary conditions Eqs. (2)and(3). Forthispurpose,domain isdividedinto nﬁnite elements, such that =nÍ e=1 e: (4) Each element ehas a relative density eand an angle e, which corresponds to the ﬁbre orientation. Thus, element e has an elasticity matrix Ce=Cee;eparametrised by Ce=P eCg e; (5) whereCg e=Cg eeistheelasticitymatrixofthebasema- terial in a global coordinate system and Pis a penalisation factor. Aplane-stressconsiderationisadopted,andtheelas- ticity tensor is Cl=b f fdC1111C11220 C1122C22220 0 0 C1212c g ge(6) in its local coordinate system and this is transformed to the global coordinate system with Cg e=RT eClRe; (7) whereRe=Reeis the rotation matrix for element e.The optimisation problem is stated as min ;f=FTU s.t.VfV KU=F minfef1:0(8) wheref=f.;/is the compliance of the structure, Fis the vector of external forces, U=U.;/is the displace- mentvector, V=V./isthevolumeofthestructure, Visa ﬁxedvolume(upperbound), K=K.;/istheglobalstiﬀ- ness matrix and minis the minimum value adopted for the relative density, used to avoid numerical issues during both the optimisation procedure and solving the linear system. TheoptimisationproblemstatedinEq. (8)canbewritten asanunconstrainedproblemusingtheLagrangianfunction, leading to min ;L=FTU+V*V s.t.KU=F minfef1:0;(9) whereL=L.;/is the Lagrangian function and is the Kuhn-Tucker multiplier. The minimum is reached when (;L=0: (10) Note that the volume constraint is dependent only on  (andindependentof ). Therefore,theprocedurecanbesep- aratedintoaconstrainedandanunconstrainedproblemtobe solved sequentially. Derivating Eq. (10) with respect to , returns (L=(f+(V=0: (11) Asshownin[30],the q-threlativedensitycanbeupdated by k+1 q=k q  q; (12) whereis a relaxation parameter and qis the update pa- rameter, given by q=df dq dV dq: (13) The derivatives of the compliance and of the volume in relation to the q*th relative density are given by df dq=uT qHpgÉ m=1BT q PP*1 q Cg qBqWmJmI uq(14) JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 3 of 15 Concurrent topology and ﬁbre angle optimisation of composites and dV dq=Vq; (15) respectively, where Bqis the strain-displacement matrix of theq-th element of the mesh, pgis the number of Gauss points used for the numerical integration, WmandJmare, respectively, the quadrature weight and determinant of the Jacobian matrix associated to Gauss point m. Now, taking the gradient from Eq.(10) with respect to , gives (L=0=(f; (16) which is an unconstrained problem. The derivative of the compliance with respect to the q-th angle is df dq=uT qHpgÉ m=1BT qP qdCg q dqBqWmJmI uq(17) in which the derivative of the elasticity matrix with respect toqcan be obtained analytically from Eq.(7). The un- constrained problem is solved using the Steepest Descent method with a Golden section algorithm for the line search. Considering that the derivatives represent a local be- haviour of the function, a scheme of moving limits is adopted. These constraints are k+1 qË k q*d1 q;k q+d1 q (18) for the relative density and k+1 qË k q*d2 q;k q+d2 q ; (19) for the angle of each element, where d1 qandd2 qare positive moving limits for the q-th relative density, and the q-th an- gle,respectively. Themovinglimitschangeateachiteration depending on the behaviour of the design variable. If the design variable changes monotonically in three subsequent steps, the moving limit is updated by a factor higher than 1:0, and if the design variable oscillates in those steps, the moving limit is updated by a factor lower than 1:0. Intopologyoptimisationproblems,twosigniﬁcantissues may occur. The ﬁrst one is the appearance of patterns simi- lartoacheckerboard,inwhicharegionhas,alternately,solid andvoidelements. Thesecondoneisthemesh-dependency oftheresults,inwhichdiﬀerentresultsareobtainedfordif- ferentmeshsizes[31]. Thus,aﬁlteringschemeisadoptedto avoid numerical instabilities in the optimisation procedure. Abasicsensitivityﬁlteringscheme[31]isusedforbothrel- ative density and ﬁbre angle. Consider a region of radius R around an element e, which is given by Ne= j;ñcj*ceñfR ; (20)wherecjandcearethecentroidofelements jande,respec- tively. The dependency of the design variable of element e on its neighbours is written as )f )e=³ jËNewcjj)f )j e³ jËNewcjvj; (21) wherew.cj/isaweightingfunction,chosenhereasalinear decaying one. Aﬂowchartofthewholeoptimisationprocedureisshown in Figure 2. All steps are performed with an in-house code written in Julia Language [32]. The data needed for the op- timisation is an input to the algorithm. This includes the geometry,mesh,allthevariablesforthetopologyoptimisa- tion, material properties and boundary conditions. The last step before starting the optimisation procedure is to create twovectors1and1containingtheinitialestimateforboth design variables. The input data used in all simulations is shown in Appendix A. The optimisation procedure itself is performed sequen- tially. For each global optimisation step, the topology op- timisation is performed ﬁrst, followed by the material opti- misation. During the topology optimisation, the equilibrium prob- lem for the current step nis solved for nandn. Then, the derivativesofthecomplianceaboutthe nareobtainedwith Eq.(14)andusedtoupdatetherelativedensitiesto n+1ac- cording to Eq. (12) while respecting the moving limits im- posed in Eq.(18). Thematerialoptimisationstepisconductedwiththevari- ablesn+1andn, see Figure 2. First, the derivatives of the compliance about the ﬁbre angles are evaluated with Eq.(17). Second, the ﬁbre angles are updated to n+1us- ingthesteepestdescentalgorithm,whilerespectingthemov- ing limits imposed in Eq.(19). In both topology and mate- rial optimisation, a sensitivity ﬁltering scheme, Eq.(21), is adopted. Convergence is reached when the compliance, relative densities, and ﬁbre angles have a variation inferior to 1~in tensubsequentsteps. Ifconvergenceisnotsatisﬁed,oneas- sumesn=n+1, and the procedure returns to the topology optimisation. Whenconvergenceisreached,aplaintextﬁle is generated containing the compliance, a vector of relative densities and a vector of ﬁbre angles. This format is such that the optimum geometry and ﬁbre orientation can be vi- sualised using the software Gmsh [33]. 3. The optimisation 3.1. Study cases The optimisation framework is tested on three design cases: cantilever, 3-point bending and 4-point bending beams,asshowninFigure3. Theseareclassicalbenchmark cases in TO, allowing us to compare our results with other JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 4 of 15 Concurrent topology and ﬁbre angle optimisation of composites Data Reading Mesh GenerationStart Define n = 1; ρ1, θ1 Conver gence?No Yes Return f, ρn+1, θn+1 Solve KU = F K = K( ρn, θn) Find df/d ρn Gradient Filtering Optimum Criterion: Find ρn+1 n = n+1 Solve KU = F K = K( ρn+1, θn) Find df/d θn Gradient Filtering Steepest Descent: Find θn+1 TOPOLOGY OPTIMISA TIONMATERIAL OPTIMISA TIONGLOBAL OPTIMISA TION Figure 2: Flowchart of the optimisation framework. approaches available in the literature. In all cases, the prob- lem is 2D (plane stress) and the optimisation domain has a heighta= 30mm and a width b= 60mm, see Figure 3. Considering that the symmetry line is located at b= 0, the supporting pins are located at 0:9bfor both 3-point and 4- pointbending,whereastheforceonthe4-pointbendingcase is applied at 0:1b. The domain is discretised by 18090 regular bi-linear isoparametric elements with incompatible modes [34]. An example of the ﬁnite element mesh, along- side the ﬁnal optimisation topology, is shown in Figure 4. A ﬁlter radius R= 0:5mm is used in all cases to guaran- teethattwoadjacentelementsarewithintheﬁlteringradius. Foreachdesigncase,theoptimisationisconductedforthree diﬀerent ﬁnal volumes: V=30~, 40%, and 50% of the ini- tial design domain. The optimisation is conducted with two diﬀerent materi- als, one isotropic and one orthotropic, to quantify the eﬀect ofﬁbreorientation. Theisotropicsolidhasthepropertiesof nylon, whereas the orthotropic material, has those of onyx. These two materials are the same as those used in the ex-Table 1 Experimentally measured elastic and physical properties for both nylon and onyx materials. Nylon (n) Onyx (o) E(GPa) 0.98 E1(GPa) 1.26 (–) 0.42 E2(GPa) 0.33 (g/cm3) 1.1G12(GPa) 0.37 G(GPa) 0.35 12(–) 0.39 (g/cm3) 1.18 periments presented in Section 4.2. The properties of both materials are presented in Table 1. The isotropic material, nylon, is characterised by an elastic modulus E=0:98GPa and a Poisson’s ratio =0:42. Ontheotherhand,Onyxisacompositematerialinwhich nylonisreinforcedwithchoppedcarbonﬁbres. Onyxcomes asaﬁlament(tobeusedinfusedﬁlamentfabrication)where the chopped ﬁbres are predominantly aligned with the ﬁla- mentdirection,whichleadstoanorthotropicbehaviour. The JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 5 of 15 Concurrent topology and ﬁbre angle optimisation of composites (a) (c) (b) Figure 3: Design domain, loads, and boundary conditions for (a) cantilever, (b) 3-point bending, and (c) 4-point bending cases. Figure 4: The convergence–free ﬁnite element mesh used in all optimisation cases. Onyxﬁlamentshaveaﬁbrevolumefraction Vfù15~,and the carbon ﬁbres have an average diameter f= 6µm and lengthlf= 100 µm. UniaxialtestshavebeenperformedinbothOnyxandNy- lonmaterialstocharacterisetheirelasticmoduli, E1foronyx andEfornylon,andthePoisson’sratio, ,fornylon. Those values are given in Table 1. With these characteristics, it is possibletoderivetheotherpropertiesfortheorthotropicma- terialusingtheHalpin–Tsaimodelforshortﬁbre-reinforced composites[35]. Thetransverseelasticmodulus E2isgiven by: E2=Em01+2Vf 1*2Vf1 ;2=E1;f Em *1 E1;f Em +;=2(22) where the elastic moduli of the ﬁbres E1;f= 230GPa and ofthematrix Em=0.98GPa. Otherwise,thein-planeshear modulusG12is given by: G12=Gm01+3Vf 1*3Vf1 ;3=G12;f Gm *1 G12;f Gm +2;=1(23) where the shear moduli of the ﬁbres G12;f= 104.5GPa (G12;f=E1;f_2.1+12;f)andofthematrix Gm=0.35GPa (Gm=Em_2.1+m). Finally, the major Poisson’s ratio isTable 2 Nomenclature for the specimens and their parameters. Case V () R (mm) Material cb_30_n 30% 0.5Nylon cb_30_o Onyx 3pb_30_n Nylon 3pb_30_o Onyx 4pb_30_n Nylon 4pb_30_o Onyx cb_40_n 40% 0.5Nylon cb_40_o Onyx 3pb_40_n Nylon 3pb_40_o Onyx 4pb_40_n Nylon 4pb_40_o Onyx cb_50_n 50% 0.5Nylon cb_50_o Onyx 3pb_50_n Nylon 3pb_50_o Onyx 4pb_50_n Nylon 4pb_50_o Onyx given by: 12=12;fVf+mVm;12=21E2 E1(24) wherethePoisson’sratiooftheﬁbres 12;f=0:1andthatof the matrix is m=0:42. The orthotropic properties of onyx are included in Table 1 whereas an overview of all design cases considered is provided in Table 2. 3.2. Optimisation results Aconvergenceanalysisoftheobjectivefunctionisshown inFigure5forthe3pb_30_noptimisationcase. Inthiscase, the compliance is calculated at the end of each global step of the optimisation procedure, meaning that both topology andmaterialoptimisationsareperformedtoobtainthevalue of the objective function. All optimisation cases follow the same convergence behaviour, therefore only one plot is pre- sented as it is representative of all cases. Asstatedearlier,convergenceisreachedwhentheobjec- tive function and the design variables do not vary by more than1~in ten subsequent global steps. Furthermore, the JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 6 of 15 Concurrent topology and ﬁbre angle optimisation of composites 02 55 07 51 001 2502x10-24x10-26x10-28x10-2Compliance (m/N)I teration Figure 5: Convergence plot for the 3pb_30_n case. convergence is monotonic, which shows the consistency of the implemented algorithm. Finally, for both isotropic and orthotropiccases,theobjectivefunctionbehavesakintothe convergence analysis shown in Figure 5. The optimal topologies for isotropic and orthotropic ma- terials are compared in Figure 6. The results show that the topologies for the isotropic and orthotropic materials diﬀer fromeachother;therearediﬀerencesinthereinforcements, especially for V= 0:3. The nature of the materials them- selvesmakesitreasonabletoexpectsuchdiﬀerences. Inthe orthotropic case, the ﬁbres have a preferential stiﬀness di- rection, and the direction of the ﬁbres is optimised along- side the topology. In this case, the topology optimisation stepinﬂuencesthematerialoptimisationstepandvice-versa. Those features are not present in the isotropic case. More- over, when the same case with diﬀerent volume fractions is considered,onecannoticeapatterninthetopologies. Basi- cally, the higher the admissible volume of the structure, the bigger the reinforcements and the larger the number of re- inforcements. Thus,theimplementedalgorithmsshowcon- sistency. In addition, Figure 7 shows the optimal ﬁbre distribution for the orthotropic cases. The ﬁbre angles follow the rein- forcements, which can be seen in details in the highlighted areas. Thiswasexpectedsincetheﬁbresorientedinthespan direction provide a higher stiﬀness to the structure, thus re- ducing its compliance. 4. Experimental details The optimal 3-point and 4-point bending beams were manufactured and tested to verify the eﬃciency of the op- timisation algorithm. The cantilever beam was excluded fromthesetestsasitisdiﬃculttoachieveperfectlyclamped boundaryconditionsinexperiments. Thissectiondetailsthe approachfollowedtoconverttheTOoutputintoanSTLﬁle, as well as the 3D printing and testing procedures.4.1. The STL ﬁle generation Mostrapidprototypingtechniquesand3Dprintersrequire the geometry to be provided as an STL ﬁle [36]. There- fore, we have to convert the output from TO into a three- dimensional STL ﬁle. This is done using three diﬀerent toolsasdepictedontheﬂowchartinFigure8. First,arough STLﬁleiscreatedusinganin-housealgorithm. Second,the STL mesh is veriﬁed and smoothed using MeshLab. Third, Blenderisusedtocorrectminorissueswiththemesh,ifnec- essary. The entire procedure is detailed below. The ﬁrst step is performed with an in-house algorithm written in Julia Language [32]. The input for the algorithm is a plain text containing the connectivity of each element of the FE mesh and its relative density. The optimisation was done on a 2D plane stress model and using symmet- ric boundary conditions (see Section 3); therefore, our al- gorithm uses extrusion and mirroring to generate the full 3D geometry. The optimisation process often produces el- ements with intermediate relative densities ( 0< e<1), which cannot be part of the ﬁnal geometry. Our algorithm uses a threshold value tsuch that elements with relative densitiesabovethethresholdareassumedassolidandthose belowtare considered voids (see Figure 1). In all cases, tisselectedsuchthattheadmissiblevolumeismaintained. Afterwards, the algorithm generates a crude STL ﬁle from optimisationdata. ForeachfaceofeveryelementontheFE mesh,atriangleofcoordinatesisgeneratedalongsideanor- malvector. Alltrianglesfollowtheright-handruleandtheir normals point outward. Next,weusethemeshprocessingsoftwareMeshLab[37] to verify and improve the STL model. The STL ﬁle created in the ﬁrst step contains duplicate faces and vertices for ad- jacent elements, and we use MeshLab to remove these du- plicate features. Then, the mesh quality is veriﬁed to en- sure that it is ’watertight’, free of holes/gaps, and does not contain any intersecting/overlapping triangles [38]. If any problemsaredetected,thenthemeshiscorrectedmanually, and this is done is two steps. First, the mesh is simpliﬁed to a two-dimensional geometry by deleting elements in the out-of-planedirection. Second,the2Dmeshisexportedand openedwiththesoftwareBlender[39],whichincludesbuilt- in tools for repairing STL meshes. The main issue encoun- tered inthis workconsisted ofelements connectedby a sin- glevertexasshowninFigure8. Aftercorrectingtheissues, the geometry is extruded into a 3D part and re-opened in MeshLab to verify the mesh quality. Once the quality of the mesh is satisfactory, the last step is to smooth the geometry to eliminate the pixelated con- tours caused by the optimisation procedure. This is done with the Laplacian smoothing method, where the position ofeachvertexisadjustedbasedontheweightedpositionsits neighbours [40]. Finally, a new STL ﬁle is generated con- tainingthequality-checked,smoothedmesh. Then,theﬁbre orientation of each element is translated to the 3D printer afterwards. JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 7 of 15 Concurrent topology and ﬁbre angle optimisation of composites cb_30_n cb_30_o cb_40_n cb_40_o cb_50_n cb_50_o 3pb_30_n 3pb_30_o 3pb_40_n 3pb_40_o 3pb_50_n 3pb_50_o 4pb_30_n 4pb_30_o 4pb_40_n 4pb_40_o 4pb_50_n 4pb_50_o (a) (b) (c) Figure 6: Optimal topologies for (a) cantilever, (b) 3-point bending, and (c) 4-point bending beams. 3pb_30_o 3pb_40_o 3pb_50_o4pb_30_o 4pb_40_o 4pb_50_o Figure 7: Optimal ﬁbre angle for the optimised (a) 3-point bending and (b) 4-point bending beams. 4.2. Additive manufacturing & testing All samples were printed using a Mark Two 3D printer (Markforged Inc., USA), which uses fused ﬁlament fabrica- tion. Both printing materials (nylon and onyx) were stored inadrystorageboxtolimitmoistureretentionpriortoprint- ing. Theﬁlamentswereheatedwithintheprinter’sheadand laid layer-by-layer, consolidating under atmospheric condi- tions. From each simulation, the .stl ﬁles were sliced using the Markforged cloud-based software, Eiger [41]. A layer height of 0.1 mm and solid inﬁll were selected to provideas accurate detail as possible. All parts had the same over- alldimensions: 120308 mm3,andweremadeupof80 layers. All samples were tested in 3-point and 4-point bending. The tests were carried out with a displacement rate of 2 mm/min in a MTS universal testing machine equipped with a load cell of 30kN. A span-to-thickness ratio of 16:1 was used,basedontherecommendationsofASTMD7264stan- dard. 4.3. Experimental results The experimental load versus displacement curves for all nylonandonyxtopologiesunder3-pointbendingareshown in Figure 9 (graphs with diﬀerent axis limits for better vi- sualisation of results and curve shapes). All samples have a fairly ductile response, which is beneﬁcial for composite structuralcomponentssincetheirinherentbrittlenessisper- haps the main bottleneck that prevents their use in primary structures. As expected, the stiﬀness and strength increase with increasing volume. This holds true for both materials. Another remarkable characteristic is the excellent repeata- bility: eachtestwasrepeatedﬁvetimesandthestandardde- viation is very low for both stiﬀness and strength. Focusing on the eﬀect of the parent material, we observe the follow- ing: 3-point bending : •Nylon (Figure 9(a,c,e)): regardless of the volume, the curves have a similar shape, with a linear–elastic in- crease up to a peak force, followed by a gradual soft- ening due to plastic deformation. Increasing the ﬁnal volume increases the deﬂection at peak force, which is mainly due to an increase in structural stiﬀness and strength. •Onyx (Figure 9(b,d,f)): in general, the responses are similar, but diﬀerent from those of nylon sam- JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 8 of 15 Concurrent topology and ﬁbre angle optimisation of composites Figure 8: Flowchart of the procedure to transform the output from optimisation into a printable STL ﬁle. ples. Firstly, these samples are signiﬁcantly stiﬀer and strongerthanthenylonones,whichwasexpectedsince theyhavereinforcingﬁbres. Nevertheless,eventhough thereinforcingﬁbresareassmallas100 µm,theyhave a strong eﬀect on the structural response of the parts thanks to the great eﬃciency of the concurrent opti- misation framework herein developed. Secondly, onyx samples exhibit a sharp load drop after the peak force, followed later by a gradual softening. No full loss of structural response is observed for any parts, which is attractive for structural components since the part can stillcarryloadaftercrackingand/orelasticbucklingoc- cur. 4-point bending : •Nylon (Figure 10(a,c,e)): all responses are similar to thosemeasuredunder3-pointbending;thereisalinear- elasticregimeuptoapeakload. Nonetheless,thepost- peakplasticdeformationisslightlydiﬀerentandoccurs ataroughlyconstantloadlevel. Thisisattributedtothe betterloaddistributionontothecompressivesideofthe sample.•Onyx (Figure 10(b,d,f)): again, these specimens are signiﬁcantly stiﬀer and stronger than nylon samples. Similarly to nylon specimens, onyx samples display very little or no post-peak softening under 4-point bending. 5. Discussion The ﬁnal deformed shapes, after unloading, for all parts are shown in Figure 11. There are no fractured trusses ob- served on any specimens. For nylon samples, plastic de- formation and local buckling are the most dominant failure mechanisms. Otherwise, minor cracking and elastic buck- ling are the main failure mechanisms for onyx specimens. Onyx parts also have more pronounced out-of-plane defor- mationandintralaminarfractureonthecompressivesideof the specimens (upper edge) when compared to nylon sam- ples. Itisworthmentioningthatinterlaminarfailure(delam- ination) was not observed in any specimens. We anticipate that delamination would be more prevalent if the polymer was reinforced with continuous ﬁbres, in-line with the ob- servations of Chen and Ye [20]. JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 9 of 15 Concurrent topology and ﬁbre angle optimisation of composites (a) (b) (c) (d) (e) (f) Figure 9: Experimental load-deﬂection curves for 3-point bending beams: (a) 3pb_30_n, (b) 3pb_30_o, (c) 3pb_40_n, (d) 3pb_40_o, (e) 3pb_50_n, (f) 3pb_50_o. A representative photograph of the part taken at the peak load is shown for each group.Graphs with diﬀerent axes limits for better visualisation of the curves. JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 10 of 15 Concurrent topology and ﬁbre angle optimisation of composites (a) (b) (c) (d) (e) (f) Figure 10: Experimental load-deﬂection curves for 4-point bending beams: (a) 4pb_30_n, (b) 4pb_30_o, (c) 4pb_40_n, (d) 4pb_40_o, (e) 4pb_50_n, (f) 4pb_50_o. A representative photograph of the part taken at the peak load is shown for each group.Graphs with diﬀerent axes limits for better visualisation of the curves. JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 11 of 15 Concurrent topology and ﬁbre angle optimisation of composites 5 5 5 5 Figure 11: Deformed shapes for all 3D printed parts after unloading. The structural stiﬀness Kis given in Figures 12(a) and 12(d)for3-pointand4-pointbendingtests,respectively. As expected, increasing the admissible volume increases the structural stiﬀness. In addition, the structural stiﬀness of Onyx specimens is considerably higher than that of nylon samples. This is surprising considering that the properties ofnylonarebetweentheboundsoftheonyxproperties,i.e., E2< E < E1, see Table 1. This improvement in struc- turalstiﬀnessshowstheeﬃciencyandpurposeofthepresent framework. The density-normalised stiﬀness is shown in Figure 12(b,e), whereas the volume-normalised stiﬀness is given in Figure 12(c,f). Both the density- and volume-normalisedstiﬀnessesconﬁrmthatonyxsamplesaresigniﬁcantlystiﬀer than their nylon counterparts. 6. Conclusions In this work, a simultaneous topology and ﬁbre orienta- tion optimisation framework has been successfully devel- oped and applied to optimised parts with isotropic and or- thotropic material properties. The objective function for all optimisation cases was to minimise compliance with three distinct volume constraints: 30%, 40%, and 50% of an ini- tialrectangulardomain. Threeclassicalbenchmarkcasesare considered: cantilever beam, 3-point bending, and 4-point JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 12 of 15 Concurrent topology and ﬁbre angle optimisation of composites Structural stiffness (a) (d) (b) (e) (c) (f)Structural stiffness Figure 12: Structural stiﬀness for specimens under (a) 3-point and (b) 4-point bending. Density-normalised stiﬀness for (c) 3-point and (d) 4-point bending tests. Volume-normalised stiﬀness for samples under (e) 3-point and (f) 4-point bending. bending. The optimised parts under 3-point and 4-point bending were 3D printed using FFF technique and tested to validatetheproposedframework. Samplesoptimisedforthe cantilever beam loading case were not considered in the ex-perimentalcampaignbecausetheboundaryconditionsused intheoptimisationscouldnotbereproducedexperimentally with the equipment available. Theproposedframeworkwasextremelyeﬀectiveatmax- JHS Almeida Jr, B Christoﬀ, V Tita, L St-Pierre: Preprint submitted to Elsevier Page 13 of 15 Concurrent topology and ﬁbre angle optimisation of composites imisingthestructuralstiﬀnessforbothnylon(isotropic)and onyx (orthotropic) parts. Experimental results showed that thedensity-andvolume-normalisedstiﬀnessesofonyxparts weresigniﬁcantlystiﬀerthanthoseofnylonsamples,which indicates that the concurrent framework was extremely ef- ﬁcient to optimise the chopped CFRP composite parts at diﬀerent admissible volumes. This framework can play a keyroleinsavingweightforadditivelymanufacturedCFRP structures. In the future, we plan to extend this framework to minimise weight under design load constraints. Acknowledgements HA is supported by the Royal Academy of Engineer- ing under the Research Fellowship scheme [Grant No. RF/201920/19/150]. VT acknowledges CNPq (process No. 310656/2018-4) and CAPES-FCT (No. AUXPE 88881.467834/2019-01 - Finance Code 001).