A concurrent fibre orientation and topology optimisation framework for 3D-printed fibre-reinforced composites

This work proposes a novel framework able to optimise both topology and fibre angle concomitantly to minimise the compliance of a structure. Two different materials are considered, one with isotropic properties (nylon) and another one with orthotropic properties (onyx, which is nylon reinforced with chopped carbon fibres). The framework optimises, in the same particular sub-step, first the topology, and second, the fibre angle at every element throughout the domain. For the isotropic material, only topology optimisation takes place, whereas, for the orthotropic solid, both topology and fibre orientation are considered. The objective function is to minimise compliance, and this is done for three volume fractions of material inside the design domain: 30%, 40%, and 50%. Two classical benchmark cases are considered: 3-point and 4-point bending loading cases. The optimum topologies are further treated and manufactured using the fused filament fabrication (FFF) 3D printing method. Key results reveal that the absolute stiffness, density-normalised and volume-normalised stiffness values within each admissible volume are higher for onyx than for nylon, which proves the efficiency of the proposed concurrent optimisation framework. Moreover, although the objective function was to minimise compliance, it was also effective to improve the strength of all parts. The excellent quality and geometric tolerance of the 3D-printed parts are also worth mentioning.


Introduction
A fundamental engineering challenge is how to design a structure to be as light as possible without sacrificing its mechanical performance.According to Sigmund and Maute [1], this can be achieved with topology optimisation (TO).Its basic idea relies on repeated analysis and design update steps usually guided by a gradient-based computation.The first attempt on TO was carried out by Bendsøe and Kikuchi [2] with the aim to propose an alternative method to shape optimisation approachable to yield both 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 function as a continuous design variable within the optimisation problem.In fact, Bendsøe and Kikuchi [2], Bendsøe [3] and Rozvany et al. [4] proposed the most disseminated mathematical approach for TO, the well-known Solid Isotropic Material with Penalisation (SIMP) method [5].This method finds the optimal material distribution within a particular design domain, load cases, boundary conditions, manufacturing constraints, and performance requirements [6,7].Other less explored approaches include, for instance, the level set [8], evolutionary structural optimisation [9], and moving morphable component methods [10].
In the last 30 years, several TO approaches have been explored, mostly for isotropic materials [11].This is mainly due to two reasons: (i) the high maturity level of approaches for isotropic materials (such as SIMP) and (ii) well-established manufacturing techniques for complex metallic shapes.However, the continuous demand for lightweight structures has led to the increased use of anisotropic carbon fibrereinforced polymer (CFRP) composites, mainly 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 filament fabrication (FFF) is one of the most disseminated AM techniques, in which a polymer filament, either reinforced with fibres or not, is extruded layer-by-layer and deposited on a build plate [15][16][17].According to the systematic review of Sanei and Popescu [18], https://doi.org/10.1016/j.compscitech.2022.109872Received 20 June 2022; Received in revised form 7 December 2022; Accepted 8 December 2022 most studies on FFF for CFRP parts were published in the last five years [19].For instance, Chen and Ye [20,21] 3D printed composite parts with carbon fibre 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 fibres placed along the principal stress directions.
However, there is still an open question: how to arrange the reinforcing fibres and optimise the part topology?Aiming at taking full advantage of both material anisotropy and exploiting the powerful manufacturing capability of FFF method, an adequate answer is by combining TO to optimise the material distribution with structural optimisation to adjust the local fibre orientation.On the one hand, the material distribution can be optimised using SIMP; on the other hand, the fibre angle can be locally optimised with gradient-based algorithms, in which the fibre angle can be treated as a design variable [23].
There are a few studies dealing with both topology and fibre angle optimisations simultaneously.For example, Lee et al. [24] developed a TO method for optimising both material layout and fibre orientation in functionally graded fibre-reinforced composites, in which the fibre angle is considered as a discrete variable.Tong et al. [25] built a sequential optimisation method considering both fibre angle and topology for constant-stiffness laminated plates using lamination parameters as design variables.The stiffness for a short cantilever beam and the flexibility for a compliant inverter increased by 6.5% and 4.2%, respectively.Jiang et al. [26] proposed a TO approach for continuous fibre angle optimisation, which computes the best layout and orientation of fibre reinforcement for AM structures.They report minimum compliance 63% lower than the baseline by selecting a different print orientation, in which the fibre 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 optimisation results have some areas which violate the conditions for tensor invariants.Papapetrou et al. [28] developed an optimisation framework for both topology and fibre paths to create variable-stiffness designs.In general, the optimised part is stiffer than the baseline.Yan et al. [29] proposed a concurrent hierarchical optimisation methodology considering the simultaneous optimisation of structural topology and orthotropic material orientation.They showed by means of numerical examples that optimising both topology and fibre angle might decrease compliance.Among these studies, only Jiang et al. [26] manufactured optimum beams using carbon-fibre-reinforced polylactic acid (PLA) composites.However, only one type of specimen and one volume were 3D-printed.Nevertheless, all the other mentioned studies are merely computational.Chen and Ye [20] developed a procedure combining topological design and fibre placement paths based on average load transmission trajectories for 3D-printed CFRP parts using the classical SIMP method for optimising the topology.
After the non-extensive state-of-the-art on research undertaken on the topic, the following gaps have been identified, which underlie the realisation of this study: • there is a need for a framework that can simultaneously optimise both topology and fibre angle; • no investigations are considering chopped carbon fibre reinforced thermoplastic composites; and • there is a lack of comprehensive investigations considering the whole design process: from optimisation to manufacturing and testing.
In this context, this work proposes a robust concurrent topology and fibre angle optimisation framework for 3D-printed composites using dedicated algorithms considering the orthotropy of the composite materials in contrast to the classical SIMP approach.The topology problem is solved as a constrained problem whereas the fibre angle is treated as an unconstrained optimisation problem.The design domain is fixed and the optimisation is done for three volume fractions of material: 30%, 40%, and 50%.In all cases, the objective function is to minimise compliance.Two benchmark cases are considered to evaluate the novel approach herein developed: (i) three-point bending (3PB), and (ii) fourpoint bending (4PB).Two material properties are considered: isotropic and orthotropic.The isotropic material consists of nylon (Polyamide 6 -nylon), whereas the orthotropic material is chopped carbon-fibrefilled nylon, called onyx.After the numerical optimisation is finished, 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 optimisation.

The formulation
Let  be a domain in R 2 with boundary  (Fig. 1).Consider that the Dirichlet and Neumann boundary conditions are applied on   and   , respectively, where such that where   is the −th component of the displacement and    is known, and where  is the stress tensor,  is a normal unit vector and   is the boundary traction vector.
We want to find the optimal material distribution that minimises compliance when subjected to boundary conditions Eqs. ( 2) and (3).For this purpose, domain  is divided into  finite elements, such that The material properties are considered constant within each element, and the variables are its relative density   and its fibre orientation angle   .Thus, an element  has an elasticity matrix   =   (   ,   ) .The effect of the relative density on the material properties is represented by a power-law [11] according to where is the elasticity matrix of the base material in a global coordinate system and  is a penalisation factor.We used  = 3 such that intermediate relative densities become unfavourable in the optimisation procedure.A plane-stress state is adopted, in which the elasticity tensor is in its local coordinate system and this is transformed into the global coordinate system with where is the rotation matrix for element .A compliance minimisation problem subject to a volume constraint is addressed here.As well established in classical works dealing with TO problems [2,3,5,11,[30][31][32], when a compliance optimisation is considered, the objective function can be calculated as the integral over a local function dependent on the displacement field and the design variables, e.g., the strain energy density.The problem consists of finding the material distribution and the angle of the fibres within the domain, such that force times displacement is minimised.The compliance minimisation problem [30,33] can thus be mathematically stated as where  =  (, ) is herein considered a compliance measure of the structure,  and  =  (, ) are the global external force and displacement vectors, respectively,  =  (, ) is the global stiffness matrix,  =  () is the volume of the structure, V is a fixed admissible volume (upper bound), and   is the minimum value adopted for the relative density, used to avoid numerical issues during both the optimisation procedure and solving the linear system.The volume constraint represents a limit on the amount of material used in the optimisation procedure.The limit is defined as V =  0   , where  0 is the volume of the design domain  and   is a prescribed allowable volume.The volume of the structure is constrained only by the relative densities.Thus, the problem dealing with the fibre angle optimisation is unconstrained, whereas the topology optimisation is constrained to the allowable volume.
The optimisation problem stated in Eq. ( 8) can be written as an unconstrained problem using the Lagrangian function, leading to where  =  (, ) is the Lagrangian function and  is the Kuhn-Tucker multiplier.The minimum is reached when The volume constraint is dependent only on  (independent of ).Thus, the procedure can be separated into a constrained and an unconstrained problem to be solved sequentially.Differentiating Eq. ( 10) with respect to , returns The relative density associated to the −th finite element can be updated by [30] where  is a relaxation parameter and   is the update parameter, given by The derivatives of the compliance and of the volume in relation to the −th relative density are given by )   (14) and respectively, where   ,   and   are, respectively, the strain-displacement matrix, the local displacement vector and the volume of the −th element,  is the number of Gauss points used for the numerical integration,   and   are, respectively, the quadrature weight and determinant of the Jacobian matrix associated to Gauss point .Now, taking the gradient from Eq. ( 10) with respect to , gives which is an unconstrained problem.The derivative of the compliance with respect to the angle of the th finite element is )   (17) in which the derivative of the elasticity matrix with respect to   can be obtained analytically from Eq. ( 7).The unconstrained problem is solved using the Steepest Descent method with a Golden section algorithm for the line search.
Considering that the derivatives represent a local behaviour of the function, a scheme of moving limits is adopted.These constraints are for the relative density and for the angle of each element, where  1  and  2  are positive moving limits for the −th relative density, and the −th angle, respectively.The moving limits change at each iteration 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.
In TO problems, two significant issues may occur.The first one is the appearance of patterns similar to a checkerboard, in which a region has, alternately, solid and void elements.The second one is the mesh dependency of the results, in which different results are obtained for different mesh sizes [31].Thus, a filtering scheme is adopted to avoid numerical instabilities in the topology optimisation procedure.The same filtering scheme is used in the fibre angles to avoid abrupt changes of orientation during the material optimisation procedure.A basic sensitivity filtering scheme [31] is used for both relative density and fibre angle.Consider a region of radius  around an element , which is given by where   and   are the centroid of elements  and , respectively.The dependency of the design variable of element  on its neighbours is written as where (  ) is a weighting function, chosen here as a linear decaying one.
A flowchart of the whole optimisation procedure is shown in Fig. 2. All steps are performed with an in-house code written in Julia Language [34].The data needed for the optimisation is an input to the algorithm.This includes the geometry, mesh, all the variables for the topology optimisation, material properties and boundary conditions.The last step before starting the optimisation procedure is to create two vectors  1 and  1 containing the initial estimate for both design variables.The input data used in all simulations is shown in Appendix.
The optimisation procedure itself is performed sequentially.For each global optimisation step, topology optimisation is performed first, followed by material optimisation.
During the topology optimisation, the equilibrium problem for the current step  is solved for   and   .Then, the derivatives of the compliance about the   are obtained with Eq. ( 14) and used to update the relative densities to  +1 according to Eq. ( 12) while respecting the moving limits imposed in Eq. (18).
The material optimisation step is conducted with the variables  +1 and   (Fig. 2).First, the derivatives of the compliance about the fibre angles are evaluated with Eq. (17).Second, the fibre angles are updated to  +1 using the steepest descent algorithm, while respecting the moving limits imposed in Eq. (19).In both TO and material optimisation, a sensitivity filtering scheme, Eq. ( 21), is adopted.
Convergence is reached when the compliance, relative densities, and fibre angles have a variation inferior to 1% in ten subsequent steps.If convergence is not satisfied, one assumes  = +1, and the procedure returns to the topology optimisation.When convergence is reached, a plain text file is generated containing the compliance, a vector of relative densities and a vector of fibre angles.This format is such that the optimum geometry and fibre orientation can be visualised using the software Gmsh [35].

Study cases
The optimisation framework is tested on two cases: 3-point bending and 4-point bending beams, as shown in Fig. 3, given that the loads and boundary conditions from the optimisations can be experimentally reproduced.These are classical benchmark cases in TO, allowing us to compare our results with other approaches available in the literature.In all cases, the problem is 2D (plane stress state) and the optimisation domain has a height  = 30 mm and a width  = 60 mm, see Fig. 3. Considering that the symmetry line is located at  = 0, the supporting pins are located at 0.9 for both 3-point and 4-point bending, whereas the force on the 4-point bending case is applied at 0.1.The domain is discretised by 180 × 90 regular bi-linear isoparametric elements with incompatible modes [36].An example of the finite element mesh, alongside the final optimisation topology, is shown in Fig. 4. A filter radius  = 0.5 mm is used in all cases to guarantee that two adjacent elements are within the filtering radius.For each design case, the optimisation is conducted for three different final volumes: V = 30%, 40%, and 50% of the initial design domain.An overview of all design cases under investigation is provided in Table 2.
The optimisation is conducted with two materials, one isotropic and one orthotropic, to quantify the effect of fibre orientation.The isotropic solid has the properties of nylon, whereas the orthotropic material has those of onyx.These two materials are the same as those used in the experiments presented in Section 4.2.The properties of both materials      For nylon: the recommendations of the ASTM D638-22 standard for type IV specimen have been followed to determine  and  using longitudinal and transverse extensometers to measure both strains.The samples were printed with a raster angle of 45 • [37].
For onyx: the ASTMD638-22 standard has been used to determine  1 ,  2 , and  12 using longitudinal and transverse extensometers.In order to measure  1 and  12 , the raster/fibre angle of the specimens is 0 • , whereas for measuring  2 , the raster angle is 90 • .To measure  12 , samples were printed with a fibre angle of ±45 • , which is common practice for long-fibre-reinforced composites according to ASTM D3518-18 standard.Here, this adaptation has been made for dog-bone samples and all specimens failed in the gauge section.
These experimentally-measured values are provided in Table 1.It is worth noting that  2 for onyx is one-third of  for nylon because to determine  2 for onyx, the fibre angle is 90 • whilst to determine  for nylon, the raster angle is 45 • .Also, the disadvantageous orientation for onyx samples printed at 90 • in relation to nylon samples printed at 45 • relative to the loading direction may promote void-opening at a higher magnitude/rate for onyx samples [38,39].

Optimisation results
A convergence analysis of the objective function is shown in Fig. 5 for the 3pb_30_o, and 3pb_30_n optimisation cases.The specimens and their parameters are given in Table 2.For the orthotropic case, the compliance is calculated at the end of each global step of the optimisation procedure, meaning that both topology and material optimisations are performed to obtain the value of the objective function.For the isotropic case, there is no material optimisation, thus, the procedure includes only of the topology optimisation step.All optimisation cases follow the same convergence behaviour, therefore only one plot is presented as it is representative of all cases.
As stated earlier, convergence is reached when the objective function and the design variables do not vary by more than 1% in ten subsequent global steps.Furthermore, the convergence is monotonic, which shows the consistency of the implemented algorithm.Finally, for both isotropic and orthotropic cases, the objective function behaves akin to the convergence analysis shown in Fig. 5.
The optimal topologies for isotropic and orthotropic materials are compared in Fig. 6.The topologies for the isotropic and orthotropic materials differ from each other; there are differences in the reinforcements, especially for  = 30%.The nature of the materials themselves makes it reasonable to expect such differences.In the orthotropic case, the fibres have a preferential stiffness direction, and the direction of the fibres is optimised alongside the topology.In this case, the topology optimisation step influences the material optimisation step and vice-versa.Those features are not present in the isotropic case.Moreover, when the same case with different volume fractions is considered, one can notice a pattern in the topologies.Basically, the higher the admissible volume of the structure, the bigger the reinforcements and the larger the number of reinforcements.Thus, the implemented algorithms show consistency.
In addition, Fig. 7 shows the optimum fibre distribution for the orthotropic cases.The fibre angles follow the reinforcements, which can be seen in detail in the highlighted areas.This is expected since the material has a longitudinal elastic modulus considerably higher than the transversal one, and thus the fibres oriented in the span direction provide a higher stiffness to the structure, reducing its compliance.In addition, no abrupt changes in the fibre angles are noticed, due to the filtering scheme used.
It is important to point out that the material properties directly affect the magnitude and shape of the global optimum.Usually, the higher the elastic properties of the material, the stiffer the optimised structures can be.Beams under 3-and 4-point bending are more dependent on  1 , therefore, this property affects the solution process to a greater extent.

Experimental details
The optimum 3-point and 4-point bending beams were manufactured and tested to verify the efficiency of the optimisation algorithm.This section details the approach followed to convert the TO output into an STL file, as well as the 3D printing and testing procedures.

The STL file generation
Most rapid prototyping techniques and 3D printers require the geometry to be provided as an STL file [40].Therefore, we have to convert the output from TO into a three-dimensional STL file.This is done using three different tools as depicted on the flowchart in Fig. 8. First, a rough STL file is created using an in-house algorithm.Second, the STL mesh is verified and smoothed using MeshLab.Third, Blender is used to correct minor issues with the mesh, if necessary.The entire procedure is detailed below.
The first step is performed with an in-house algorithm written in Julia Language [34].The input for the algorithm is a plain text file 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 symmetric boundary conditions (see Section 3); therefore, our algorithm uses extrusion and mirroring to generate the full 3D geometry.The optimisation process often produces elements with intermediate relative densities (0 <   < 1), which cannot be part of the final geometry.Our algorithm uses a threshold value   such that elements with relative densities above the threshold are assumed as solid and those below   are considered voids (see Fig. 1).In all cases,   is selected such that the admissible volume is maintained.Afterwards, the algorithm generates a crude STL file from optimisation data.For each face of every element on the FE mesh, a triangle of coordinates is generated alongside a normal vector.All triangles follow the right-hand rule and their normals point outward.
Next, we use the mesh processing software MeshLab [41] to verify and improve the STL model.The STL file created in the first step contains duplicate faces and vertices for adjacent elements, and we use MeshLab to remove these duplicate features.Then, the mesh quality is verified to ensure that it is 'watertight', free of holes/gaps, and does not contain any intersecting/overlapping triangles [42].If any problems are detected, then the mesh is corrected manually, and this is done is two steps.First, the mesh is simplified to a two-dimensional geometry by deleting elements in the out-of-plane direction.Second, the 2D mesh is exported and opened with the software Blender, which includes builtin tools for repairing STL meshes.The main issue encountered in this work consisted of elements connected by a single vertex as shown in Fig. 8.After correcting the issues, 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 contours caused by the optimisation.This is done with the Laplacian smoothing method, where the position of each vertex is adjusted based on the weighted positions of its neighbours [43].Finally, a new STL file is generated containing the quality-checked, smoothed mesh.Then, the fibre orientation of each element is translated to the 3D printer afterwards.

Additive manufacturing & testing
All samples were printed using a Mark Two 3D printer (Markforged Inc., USA), which uses fused filament fabrication.Both printing materials (nylon and onyx) were stored in a dry storage box to limit moisture retention prior to printing.The filaments were heated within the printer's head and laid layer-by-layer, consolidating under atmospheric conditions.From each simulation, the .stlfiles are sliced using the Markforged cloud-based software, Eiger.A layer height of 0.1 mm and solid infill were selected to provide as accurate detail as possible.All parts had the same overall dimensions: 120 × 30 × 8 mm 3 , and were made up of 80 layers.The fibre angle (Fig. 7) is controlled element-wise in the STL file and translated to the slicing software.
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 an MTS universal testing machine equipped with a load cell of 30 kN.A spanto-thickness ratio of 16:1 was used, based on the recommendations of the ASTM D7264-21 standard.

Experimental results
The experimental load versus displacement curves for all nylon and onyx topologies under 3-point bending are shown in Fig. 9 (graphs with different axis limits for better visualisation of results and curve shapes).A general observation is that all samples have a fairly ductile response, which is intrinsically related to the thermoplastic matrix.For parts under 3-point bending, the structural stiffness increase was 282%, 282% and 165% for onyx parts over nylon ones for admissible volumes of 30%, 40% and 50%, respectively; whereas for parts under 4-point bending, stiffness increases of 169%, 62% and 137% were reached for onyx parts over nylon ones for admissible volumes of 30%, 40% and 50%, respectively.As expected, the stiffness and strength increase with increasing volume.This holds true for both materials.Another remarkable characteristic is the excellent repeatability: each test was repeated five times and the standard deviation is very low for both stiffness and strength.Focusing on the effect of the parent material, we observe the following:

3-point bending:
• Nylon (Fig. 9(a,c,e)): regardless of the volume, the curves have a similar shape, with a linear-elastic increase up to a peak force, followed by a gradual softening due to plastic deformation.Increasing the final volume increases the deflection at peak force, which is mainly due to an increase in structural stiffness and strength.• Onyx (Fig. 9(b,d,f)): in general, the responses are similar but different from those of nylon samples.Firstly, these samples are significantly stiffer and stronger than nylon ones, which was expected since they have reinforcing fibres.Nevertheless, even though the fibres are as small as 100 μm, they have a strong effect on the structural response of the parts thanks to the great efficiency of the concurrent optimisation 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 still carry load after cracking and/or elastic buckling occur.

4-point bending:
• Nylon (Fig. 10(a,c,e)): all responses are similar to those measured under 3-point bending; there is a linear-elastic regime up to a peak load.Nonetheless, the post-peak plastic deformation is slightly different and occurs at a roughly constant load level.This is attributed to the better load distribution onto the compressive side of the sample.
• Onyx (Fig. 10(b,d,f)): again, these specimens are significantly stiffer and stronger than nylon samples.Similarly to nylon specimens, onyx samples display very little or no post-peak softening under 4-point bending.

Discussion
The final deformed shapes, after unloading, for all parts are shown in Fig. 11.There are no fractured trusses observed on any specimens.For nylon samples, plastic deformation and local buckling are the most dominant failure mechanisms.Otherwise, minor cracking and elastic buckling are the main failure mechanisms for onyx specimens.Onyx parts also have more pronounced out-of-plane deformation and intralaminar fracture on the compressive side of the specimens (upper edge) when compared to nylon samples.It is worth mentioning that interlaminar failure (delamination) was not observed in any specimens.We anticipate that delamination would be more prevalent if the polymer was reinforced with continuous fibres, in line with the observations of Chen and Ye [20].
The structural stiffness  is given in Fig. 12(a) and (d) for 3-point and 4-point bending tests, respectively.As expected, increasing the admissible volume increases the structural stiffness.In addition, the structural stiffness of Onyx specimens is considerably higher than that of nylon samples.This is surprising considering that the properties of nylon are between the bounds of the onyx properties, i.e.,  2 <  <  1 , see Table 1.This improvement in structural stiffness shows the efficiency and purpose of the present framework.
The density-normalised stiffness is shown in Fig. 12(b,e), whereas the volume-normalised stiffness is given in Fig. 12(c,f).Both the density-and volume-normalised stiffness values confirm that onyx samples are significantly stiffer than their nylon counterparts.

Comparison between experiments and simulations
The numerical stiffness is evaluated directly through the in-house code, mimicking the procedure used in experiments.A comparison between the numerical and experimental results is shown in Table 3.A good agreement between predictions and experiments is achieved; in fact, the largest difference is only 13% observed for the 3pb_40_o case.Additionally, the experimental results slightly overestimate the numerical results for most of the cases, except for cases 3pb_30_n, 3pb_40_n, and 4pb_30_n.For all cases, the predictions are within the experimental average ±1 SD, revealing an excellent agreement between them.

Parametric analysis
A parametric study highlights the benefits of considering fibre orientation in optimisation and printing.The optimum sample for an allowable volume of 30% under 3-point bending (3pb) and 4-point bending (4pb) is selected to evaluate the following scenarios: • Scenario_1: The optimum cases 3pb_30_o and 4pb_30_o (Fig. 6) are 3D-printed with onyx material controlling the fibre angle as in Fig. 7. • Scenario_2: The optimum topologies 3pb_30_n and 4pb_30_n (Fig. 6) are 3D-printed with nylon material.• Scenario_3: The optimum cases 3pb_30_o and 4pb_30_o (Fig. 6) are 3D-printed with onyx without controlling the fibre angle.The slicing software automatically determines the nominal fibre angle of 45 • throughout the domain.
• Scenario_4: The optimum topologies 3pb_30_n and 4pb_30_n (Fig. 6) are 3D-printed with onyx material instead of nylon.As the optimum topologies for the isotropic case do not have fibre angle information, these parts are 3D-printed with onyx material with all fibres oriented at 45 • .• Scenario_5: The optimum topologies 3pb_30_o and 4pb_30_o (Fig. 6) are 3D-printed with nylon material.
As shown in Table 4, Scenario_1 has the best structural stiffness over all other scenarios.This is due to the concurrent topology and material optimisation, which takes advantage of the higher elastic modulus in a preferential direction.When the same topology is tested without controlling the fibre orientation (Scenario_3) and using nylon (Scenario_5), a significant drop in the structural stiffness is observed compared to Scenario_1.Moreover, note that Scenario_2 has a better structural response than Scenario_5.As the topology obtained for Sce-nario_2 is optimum regarding an isotropic material, it is expected that any other topologies provide an inferior structural response, which is seen in our results.In Scenario_2, and Scenario_4, the same topologies are tested using nylon material, and onyx without fibre direction control, respectively.It is expected that the topology printed with nylon presents better results, which is observed for the 4pb case.For the 3pb case, this behaviour is not observed.It might be explained by the fact that the main reinforcement of the 3pb_30_n topology is mostly oriented at 45 • , which is the nominal fibre angle throughout the domain, which is favourable to the overall stiffness of the structure.

Conclusions
In this work, a simultaneous topology and fibre orientation optimisation framework has been successfully developed and applied to optimise parts with isotropic and orthotropic material properties.The objective function for all optimisation cases was to minimise compliance (strain energy) with three distinct volume constraints: 30%, 40%, and 50% of an initial rectangular domain.Two classical benchmark cases are considered: 3-point bending, and 4-point bending.The optimised parts were 3D printed using FFF technique and tested to validate the proposed framework.
The optimisation framework was extremely effective at maximising the structural stiffness for both nylon (isotropic) and onyx (orthotropic) parts.Experimental and computational results showed that the stiffnesses of onyx parts were higher than those of nylon samples, which indicates that the concurrent framework was efficient to optimise the chopped CFRP composite parts at different admissible volumes.Considering that  11 (onyx) is 29% higher than  (nylon), it is here concluded that a framework that considers the local fibre angle optimisation is essential to enhance the structural performance of fibre-reinforced structures as an additional step to topology/shape optimisation.
An experimental parametric analysis with several scenarios has been conducted and it was concluded that considering and controlling the locally optimised fibre angle in the 3D-printed parts provides the highest possible stiffness amongst all scenarios herein carried out to test the effectiveness of the framework.Moreover, a comparison between numerical and experimental stiffnesses shows a good agreement, in which all predictions are within the experimental error (average ±1 SD).This confirms that the concurrent optimisation approach hereby developed is efficient from design to manufacturing and testing of complex 3D-printed parts.Hence, this optimisation procedure may be used to reliably address both material and fibre angle for orthotropic materials and the procedure may be expanded to tackle problems with different design domains, materials, and boundary conditions.This work represents progress in the state-of-the-art by addressing the anisotropic behaviour of carbon fibre-reinforced parts in the topology optimisation of anisotropic solids by adding a sequential sub-step to optimise the fibre angle, which is vital to increase performance whilst decreasing weight.-moving limit update percentage (lower) # deltasup -moving limit update percentage (upper) # LMinf _ theta -minimum allowed moving limit # (fibre angle) # LMsup _ theta -maximum allowed moving limit # (fibre angle) # LMinf _ rho -minimum allowed moving limit (density) # LMsup _ rho -maximum allowed moving limit (density) # stop -percentage of design variables # variation to assume convergency # stop _ iter -subsequent steps to assume convergency # const deltainf = 0.7 const deltasup = 1.2\scriptsize const LMinf _ theta = 0.01 const LMsup _ theta = 0.2 const LMinf _ rho = 0.01 const LMsup _ rho = 0.1 const stop = 0.01 const stop _ iter = 10 return simp, eta, step, tol, lambda1, lambda2, frac, N, deltainf, deltasup, LMinf _ theta, LMsup _ theta, LMinf _ rho, LMsup _ rho, Rmax, nvmax, theta _ inicial, rho _ min, rho _ max, stop, stop _ iter end #function

Fig. 1 .
Fig. 1.Schematic representation of both topology and fibre orientation within a given domain.

Fig. 4 .
Fig. 4. The convergence-free finite element mesh used in all optimisation cases.
the filament direction, which leads to an orthotropic behaviour.The Onyx filaments have a fibre volume fraction   ≈ 15%, and the carbon fibres have an average diameter   = 6 μm and length   = 100 μm.Uniaxial tensile tests have been performed in both Onyx and Nylon materials to characterise their elastic properties: longitudinal elastic modulus  1 , transverse elastic modulus  2 , in-plane shear modulus  12 , and major Poisson's ratio  12 for onyx; and  and the Poisson's ratio, , for nylon.

Fig. 8 .
Fig. 8. Flowchart of the procedure to transform the output from optimisation into a printable STL file.

Fig. 10 .
Fig. 10.Experimental load-deflection 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 different axes limits for better visualisation of the curves.

Table 1
Experimentally measured elastic and physical properties for both nylon and onyx materials.
are presented in Table1.The isotropic material, nylon, is characterised by an elastic modulus  = 0.98 GPa and a Poisson's ratio  = 0.42.On the other hand, Onyx is a composite material in which nylon is reinforced with chopped carbon fibres.Onyx comes as a filament (to be used in fused filament fabrication) where the chopped fibres are

Table 2
Nomenclature for the specimens and their parameters.

Table 3
Experimental and numerical stiffness values.

Table 4
Parametric analysis to benchmark optimised structures with both materials and with/without controlling the fibre orientation.