Spatially optimised fibre-reinforced composites with isosurface-controlled additive manufacturing constraints

A design approach accounting for manufacturing constraints is described for spatially optimised fibre-reinforced composites. The approach is based on the optimisation of local fibre orientation, the fibre volume fraction and density-based topology optimisation to determine the optimal design. A continuity equation is adopted to constrain the fibre orientation and ensure continuous fibres within the bounds of realistic fibre volume fractions. This results in a fibre orientation with a corresponding and controllable variation of the fibre volume fraction. In order to ensure the continuous fibre can be deposited, the manufacturability of the optimised results is ensured by introducing constraints controlled with two scalar fields to reconstruct fibre paths which are able to provide sufficient information to generate printer toolpaths. A cantilever beam problem is solved to show the advantage of the fibre reinforcement, the inclusion of manufacturing constraints and the penalty in compliance due to the application of the manufacturing constraints. The results show that the presented approach successfully guarantees the manufacturability with minimal loss of performance.


Introduction
Fibre-reinforced composites typically possess high stiffness, high strength and low density properties. Most fibre-reinforced composites are manufactured by wet layup, filament winding or compression moulding. Additive manufacturing has extended the potential of fibre-reinforced composites due to the automatic manufacturing process reducing the cost of labour and the precise localisation of multiple axes moving print heads increasing the ability to handle complex geometries. In additive manufacturing applications of fibrereinforced composites, fibre orientations are guided by the print head and follow deposition directions. The possibility of steering the fibre orientation spatially to fully utilise the anisotropic properties has attracted attention from researchers. As manufacturing techniques developed, early studies started to design variable stiffness laminates by optimising the spatially varying fibre orientation (Brampton et al. 2015;Khani et al. 2017). To further improve the mechanical performance and gain an extra advantage of lightweighting, the development of novel design approaches has focussed on combining topology optimisation with fibre orientation optimisation (Nomura et al. 2015;Lee et al. 2018;Boddeti et al. 2020;Caivano et al. 2020;Papapetrou et al. 2020;Schmidt et al. 2020;Jantos et al. 2020). And recently, Nomura et al. (2022) have demonstrated the optimised design not only in numerical simulations but also in practical experiments. Brampton et al. (2015) applied a level-set function to describe the equal spacing fibre path which ensures the manufacturability of laminate composites using the automated fibre placement technique. Nomura et al. (2015) adopted a topology optimisation method using isoparametric projection to generate continuous and discrete fibre orientation design. Papapetrou et al. (2020) demonstrated density-based and level-set topology optimisation of varying fibre orientation in 2D. To ensure the continuity of the fibre orientation, a fibre-filtering scheme was presented. To connect design with manufacturing, three infill methods-the equal space method, the streamline method and the offset method-were introduced to interpret the fibre path. However, the results from Brampton et al. (2015) and the offset method in Papapetrou et al. (2020) show sharp turns in the fibre path, and the results from Nomura et al. (2015) and Papapetrou et al. (2020) show varying distances between the fibre path which potentially lead to the generation of voids. The sharp turns at the fibre path result in high curvature of the deposited fibres which may cause twisting and folding back of the fibres (Matsuzaki et al. 2018) and the voids result in decreasing the fibre volume fraction. As a result, these may degrade the mechanical properties of the fibre-reinforced composites, and the neglected reduction of the fibre volume fraction may lead to a loss of modelling accuracy. Inspired by the innovation of the robotic 3D printer, Schmidt et al. (2020) presented an approach which extends the variation of the fibre orientation to three dimensional space. Their approach shows the ability to solve large-scale and multiple load optimisation problems. Despite the fact that the optimised GE bracket (GrabCAD 2013; Whalen et al. 2021) in their work achieved by the streamline method shows a smoothly varying fibre orientation which solves the problem of the sharp turns, varying distances between the fibre paths can still be observed. This illustrates that the fibre volume fraction variation corresponding to the varying distances was not considered in the modelling. In addition, it is difficult to extract sufficient information, i.e. layer height and printing width, from the three dimensional streamlines to slicing software before printing. The sharp turns and the varying distances illustrate the challenges of including manufacturability in the optimisation process.
Most applications of continuous fibres in additive manufacturing use filament deposition modelling. The fibre impregnation process directly affects the fibre volume fraction. Controlling and increasing the fibre volume fraction over 50% is a challenge (Liu et al. 2020). By pre-impregnating fibres within a matrix, the fibre volume fraction can be enhanced and the quality of impregnation can also be improved (Dickson et al. 2020). However, the fibre volume fraction is not adjustable during the deposition process. One potential method of manufacturing varying fibre volume fraction composites is using in situ impregnation. Recently, Pelzer and Hopmann (2021) have presented an algorithm for non-planar path planning with variable layer height. The variable layer height is achieved by adjusting the amount of extruded material. Adapting this concept to in situ impregnation and assuming the fibres are perfectly impregnated, as the amount of the extruded matrix varies whilst the amount of the fibre stays constant, the fibre volume fraction will have a respective variation. The assumed integration shows the manufacturability of fabricating fibre-reinforced composites with a varying fibre volume fraction. Therefore, a corresponding design approach is needed. This paper describes a novel approach to introduce manufacturing constraints to the process of designing fibrereinforced composites. The constraints prevent sharp fibre turns and generate fibre path with sufficient information for slicing software in three dimensional space. In addition, the accuracy of modelling the fibre-reinforced composites is increased by considering the fibre volume fraction as a parameter of mechanical properties. To demonstrate the efficacy of the design approach, a numerical example of a classic cantilever beam problem is presented. The topology layout, the fibre orientation and the compliance are compared for five cases to illustrate the advantage of fibre reinforcement, the interaction between the fibre orientation and the fibre volume fraction, the effectiveness of the manufacturing constraints and the loss in terms of mechanical performance caused by applying the manufacturing constraints.

Methods
The microscale model, manufacturing constraints and optimisation formulations are described in this section. The structural analysis is considered within the theory of linear elasticity, such that where ij denotes the Cauchy stress tensor, F i denotes forces acting on the design domain, E ijkl denotes the elasticity tensor, kl denotes the strain tensor and u i , u j represent the displacement.
It should be noted that the proposed method is based on a potential manufacturing method of varying fibre volume fraction composites using in situ impregnation technology which is still under development.

Microscale model
The microscale model is a unidirectional fibre-reinforced composite which can be oriented in three dimensional space. The unidirectional fibre-reinforced composite is a transversely isotropic material composed of fibre and matrix. There are different approaches to obtain the mechanical properties of unidirectional fibre-reinforced composites, such as the well-known rule of mixtures, the Eshelby equivalent homogeneous inclusion approach and the approach of representative volume elements modelled using the finite elements method (Clyne and Hull 2019). In order to keep the simplicity of the optimisation approach, the rule of mixtures using the slab model under the equal stress and equal strain assumptions is applied to evaluate the elasticity characteristics as where E denotes Young's modulus, V denotes the volume fraction, subscripts 1 and F denote the fibre direction, subscripts 2, 3 and T denote transverse directions and subscripts f, m denote fibre and matrix respectively. To capture the variation of the fibre volume fraction during optimisation, the fibre volume fraction is chosen to be one of the optimisation variables. Substituting Eq.
(2) to the elasticity tensor and rewriting it as a second order tensor gives the stiffness matrix of the microscale model The fibre orientation is defined as a vector, U = (U x , U y , U z ) . To evaluated the stiffness matrix of the microscale model with the fibre orientation U , which is not aligned with the x axis, a local coordinate system, x"y"z", in which the fibre orientation vector U can be written as (U x'' , 0, 0) , is introduced. In the x"y"z" coordinate system, the stiffness matrix can easily be evaluated by Eq. (3). By transforming the local coordinate system twice, the stiffness matrix corresponding to the global coordinate system, xyz, can be evaluated. The coordinate system transformation illustrated in Fig. 1 is achieved by rotating the local coordinate system, x"y"z", about the z" axis and then rotating the transit coordinate system, x ′ y ′ z ′ , about the y ′ axis. The rotation angles, and , can be written in terms of the components of U as The stiffness matrix after first coordinate transformation can be written as K = TK 0 T T where T is the transformation matrix (Clyne and Hull 2019). Hence, the stiffness matrix of the microscale model after second coordinate transformation can be written as a function of the fibre volume fraction and the fibre orientation vector as Note that to improve the uniqueness of the vector presented fibre orientation, the fibre orientation vector is constrained to be a unit vector, ‖U‖ = 1 , and the 2 ambiguity commonly seen in continuous fibre angle orientation optimisation (Bruyneel and Fleury 2002) is avoided by evaluating sin and cos functions directly using the component of the fibre orientation vector. However, approaches using coordinate rotations inevitably face a singular point problem. In this study, a potential singular point exists where the fibre orientation vector is aligned with the y axis, (0,1,0). This with c 1 = cos and s 1 = sin x z x y = y φ θ Fig. 1 The fibre orientation vector U is defined as (U x , U y , U z ) in the coordinate system xyz and it can be written as (U x'' , 0, 0) in the coordinate system x"y"z". The local coordinate system, x"y"z", is rotated about the z" axis through an angle which results in a transit coordinate system x ′ y ′ z ′ . The transition coordinate system is then rotated about the y ′ axis through an angle which results in the global coordinate system xyz might raise a numerical instability issue of zero division in Eq.(4), and lead to a loss of sensitivity information with respect to the second rotation at this point.

Topology optimisation
In order to gain the advantage of lightweighting and fully utilise the anisotropic properties, topology optimisation is applied in combination with fibre orientation optimisation. The most widely used density-based approach, Solid Isotropic Material with Penalisation (SIMP), is applied in this paper to describe the topology. A parameter p is used to penalise the local material stiffness matrix as where is a dimensionless pseudo density which provides a continuous interpretation of the existence of materials and min is a small number to avoid numerical instability.
Combining the penalisation with the microscale model, the pseudo density, the fibre volume fraction and the fibre orientation can be subjected to different objective optimisation problems.

Manufacturing constraints
The manufacturing problems discussed in Sect. 1 include the sharp fibre turns, the degradation of mechanical properties and the loss of modelling accuracy due to the neglected reduction of the fibre volume fraction, and the difficulty of extracting sufficient information for slicing software from streamlines. Continuity constraints and fibre path generation constraints are introduced to avoid the corresponding problems.

Continuity constraints
To overcome the manufacturing problem with regard to the fibre volume fraction and to keep the fibre orientation continuous, a continuity constraint is applied to ensure the conservation of fibre mass in an element such that where F V F denotes the fibre mass and V denotes the volume of the element. Equation (7) illustrates that the mass enters and leaves the element is the same which is analogous to the continuity equation for the case of a steady state fluid.
Removing the density of the fibre from the divergence operator since it is a constant and replacing the fraction of the fibre volume and the element volume by the fibre volume fraction gives the equation of the continuity constraint as Note that the fibre volume fraction contains the physical meaning of mass variation under the continuity constraint due to the fibre orientation vector is constrained to a unit vector.

Fibre path generation constraints
The fibre path generation constraints are inspired by fieldbased toolpath generation methods which ensure the gradient of the scalar field is orthogonal to the deposition direction to generate toolpaths (Fang et al. 2020;Chen et al. 2022). Here, the fibre path needs to be aligned with the fibre orientation vector to be consistent with the modelling result. By introducing two scalar fields S 1 and S 2 and constraining the cross product of their gradients to be equal to the fibre orientation vector as the fibre orientation vector U is orthogonal to ∇S 1 and ∇S 2 due to the orthogonality of cross products. In other words, the fibre orientation vector lies in the isosurfaces of both S 1 and S 2 . Hence, by extracting the intersection curves of the isosurfaces, the fibre path is generated. In addition, the variation of the fibre volume fraction needs to be considered in the fibre path generation. By scaling the fibre orientation vector U to V f U and rewriting Eq. (9) as the cross product of ∇S 1 and ∇S 2 is given the physical meaning of the fibre volume fraction. The magnitude form of Eq. (10) can be written as where is the angle between the gradients. To eliminate the influence of the varying angle, the gradients are constrained to be orthogonal by Since the magnitude of the fibre orientation vector is constrained to 1 and sin becomes unitary, the product of the magnitudes of the gradients is equal to the fibre volume fraction. Specifically, the fibre volume fraction is interpreted by the magnitudes of the gradients. Figure 2 shows the cross section of deposited fibre bundles where dashed lines denote the isosurfaces, subscripts 1 and 2 denote the value corresponding to S 1 or S 2 respectively, and ′ denote the scalar value of the isosurface and D denotes the distance between the isosurfaces. Assuming the voids are negligible, the fibre volume fraction can be derived as where n denotes the number of fibres within a deposited bundle and r f denotes the radius of a single strand of fibre. To extract the isosurfaces evenly, the segment values Δ used to extract isosurfaces should be the same such that Equation (14) can be written as Moving D 2 to the left hand side gives The distances between isosurfaces can be evaluated by solving Eq. (13) and Eq. (15) such that Equation (16) illustrates the aspect ratios of the deposited bundles are controlled by the ratio of the gradients. In reality, the aspect ratio of deposited fibre bundles is restricted by the nozzle geometry and the pressure in the print head. Zelinski (2022) reported an approach applying a variable width rectangular nozzle and controlling the amount of deposited matrix to decouple the link of layer height and printing width. The magnitude of the gradients can be used to control the aspect ratio in the future when the rectangular nozzle is applied to the continuous fibre filament deposition modelling. In this paper, the aspect ratio is constrained to 1 to satisfy the manufacturability requirement of current processes as

Optimisation formulations
The mechanical properties of the model are controlled by the fibre volume fraction, the fibre orientation and the pseudo density. In order to ensure the manufacturability, two scalar fields, S 1 and S 2 , the continuity constraints, Eq. (8), and the fibre path generation constraints, Eq. (10), Eq. (12) and Eq. (18), are introduced to limit the variation of the parameters of the mechanical properties. The sets of the continuity constraints C 1 and the fibre path generation constraints C 2 can be written in the form of inequality constraints as where Ω denotes the design domain. Note that all constraints are scaled by the pseudo density , which means constraints are only considered when material is present. In implementing the integrated optimisation problem, the fibre volume fraction, the fibre orientation, the pseudo density and the scalar fields are implemented as design variables. The optimisation formulations with manufacturing constraints are written as The cross section of deposited fibre bundles shows the aspect ratio of the bundles is controlled by the distances between isosurfaces, D 1 and D 2 where the objective function J is defined as the compliance, Γ denotes the domain where the force is applied, f denotes the force, u denotes the displacement, K denotes the global stiffness matrix, and V 0 denotes the target volume of the optimised result. The bound of the fibre volume fraction is set from 0 to 0.5 which models the scenario from no fibre to the average maximum value of fibre-reinforced composites manufactured by filament deposition modelling. In order to avoid potential checkerboard issues and to generate smoothly varying fibre orientations, the design variables are passed through a Helmholtz-type filter (Lazarov and Sigmund 2011) where denotes the set of the design variables, r denotes the filter radius and ̂ denotes the set of filtered design variables. To avoid the boundary effect discussed in Clausen and Andreassen (2017), the boundary condition used to solve the filtered pseudo density is changed from the homogeneous Neumann boundary condition to the Robin boundary condition which is consistent with the physical padding techniques (Wallin et al. 2020). By solving the partial differential equation, the filtered design variables are used to evaluate the rest of the optimisation problem. As the implementation of the variable filter, a chain rule corresponding to the filter is applied to evaluate sensitivities as where dJ d̂ is evaluated by the adjoint method which is implemented by the open-source Python library dolfin-adjoint (Mitusch et al. 2019). The optimisation is performed using IPOPT, a primal-dual interior-point algorithm with a filter line-search method for nonlinear programming (Wächter and Biegler 2005), and the second-order derivatives are evaluated using L-BFGS. The termination criteria are when the objective function decreases less than 0.01% and the constraint violation is less than 10 −4 in five successive iterations.

Results
In this section, a numerical example of a cantilever beam is presented to illustrate the spatially optimised fibre-reinforced composites and the manufacturing constraints. Five scenarios based on different settings for the fibre volume fraction, the fibre orientation and the manufacturing constraints are considered (see Table 1). The design variables are spatially varying within their bound constraints in the design domain except for the unidirectional case and isotropic case. For the unidirectional case, the fibre volume fraction is fixed to its upper bound and the fibre orientation is changed from spatially varying vectors to a single vector to represent a unidirectional fibre-reinforced composite. Note that the mechanical properties are evaluated by substituting V f and (23) dJ d = dJ d̂dT able 1 Scenarios for the numerical example. The design variables, V f and U , and additive manufacturing constraints are considered as control variables of the numerical example. The isotropic case uses homogeneous mechanical properties listed in Table 2 to represent a random-oriented short fibre-reinforced composite with a fixed 50% fibre volume fraction. In the unidirectional case, the fibre volume fraction is fixed to 50% and the fibre orientation is changed from spatially varying vectors to a single vector to represent a unidirectional fibre-reinforced composite. The continuity constraints are enabled in the last two cases, and the fibre path generation constraints are only enabled in the last case  (Sanomura and Kawamura 2003;Avanzini et al. 2022) to represent a random-oriented short fibre-reinforced composite with a fixed 50% fibre volume fraction. The continuity constraints are enabled in the last two cases, and the fibre path generation constraints are only enabled in the last case.

Example setup
The geometry of the cantilever beam is shown in Fig. 3. A uniform force is applied to the bottom-right as shown and the target volume of optimised results is set to the 12% volume of the design domain. A mesh convergence study concludes that discretizing the geometry into 65,887 tetrahedral elements using a mesh size of 1/30 could achieve the balance between accuracy (< 0.5% error) and the computational cost. The material properties and the initial guesses of the design variables are listed in Table 2 and Table 3. Note that several initial guesses have been tested. The convergence performance is significantly quicker when the initial guesses of the design variables are located in feasible design space. The values listed in Table 3 provide a neutral performance of the compliance and avoid large constraint violation. The penalisation factor is set to 3. The regularisation length scale of the Helmholtz filter radius is equal to the mesh size. The boundary length scale for the Robin boundary condition to solve the Helmholtz-type filter is equal to the filter radius. Note that all the present result geometry isosurfaces have threshold of 50% on for a better visualisation. The fibre orientation is visualised by streamlines with random seed points and coloured in the magnitude of the fibre volume fraction.

Cases without manufacturing constraints
The optimised topologies are shown in Fig. 4 where the isotropic case shows a typical result of topology optimisation, the topology of the unidirectional case is composed of a diagonal member and a horizontal member without holes (see Fig. 4b), and the spatially optimised case is composed of several members with three holes at the bottom-right of the topology (see Fig. 4c). The comparison of the topologies shows that parameterising the fibre orientation change the topologies. The comparison of the fibre orientation presented in Fig. 5a shows that the unidirectional case has a uniform fibre orientation pointing to the direction (9.989, 0, −1.544) and the spatially optimised case presents a varying fibre orientation parallel to the topology layout. The fibre orientation is aligned with the principal stress direction. Due to the limitation of the single orientation, the fibre orientation is only aligned with the principal stress direction of the main diagonal member in the unidirectional case. To overcome the loss of not following the principal stress direction, the horizontal member is thicker than the diagonal member (see Fig. 4b.i and b.ii).
From the manufacturability point of view, the unidirectional case can be manufactured by layering up with support materials on the xz plane. In contrast, the result shown in Fig. 5b has shown several issues that are not favourable for fabrication. First, sharp fibre turns are observed at the intersection area where the maximum principal stress direction changes. Some of the fibres connect together to form sharp turns (see Fig. 5b.i) and the rest of them result in discontinuous fibres (see Fig. 5b.ii). Second, gaps and the variation of distances between fibres can be observed (see Fig. 5b.iii) although the colour legend shows that the fibre volume fraction is almost constant. Furthermore, it is difficult to interpret the fibre volume fraction by these streamlines controlled by random seeding. Table 4 illustrates that the compliance shows significant improvements in the cases with fibres, compared to the isotropic material case. Furthermore, the spatially optimised case shows an extra 8.47% reduction, compared to the unidirectional case. Although the spatially optimised  Table 2 Material properties used in the numerical example. Note that the properties of the random-oriented short fibre-reinforced composite are used in the isotropic case Young's modulus (GPa)

Poisson ratio
Fibre 3.6 0.35 Random-oriented short fibre-reinforced composite 9.1 0.32 Table 3 Initial guesses used in the numerical example. Note that S 1 and S 2 are assigned to the coordinate value of y and z, respectively, to make ∇S 1 = (0, 1, 0) and ∇S 2 = (0, 0, 1)

Values
(1, 0, 0) 0.5 0.12 y z case shows the best performance in terms of the compliance, the sharp fibre turns may result in the degradation of mechanical properties as discussed in previous section. These may also lead to fibre breakages and therefore increase the chance of composite failure. This comparison demonstrates the advantage of the fibre reinforcement and the need for the manufacturing constraints to discourage the problematic features and bridge the gap between optimal design with manufacturing.

Cases with manufacturing constraints
The first set of the manufacturing constraints considered is the continuity constraints (see Eq. (19)). The topology shown in Fig. 6a is composed of a diagonal member and a horizontal member which is similar to the unidirectional case. The fibre orientation shown in Fig. 7a illustrates that the fibre orientation is successfully constrained by the adapted continuity equation to form a smooth transition at the intersection area without the sharp fibre turns (see Fig. 7a.i). In addition, the fibre volume fraction variation is observed at the right of the topology where the fibre orientation varies significantly. This illustrates how the fibre volume fraction is linked with the fibre orientation by the continuity constraints. As shown in Fig. 7a, the view in the xz plane shows the separation space between the fibres at the bottom right of the optimised topology decreases in the z direction whilst the view in the yz plane shows it increases in the y direction at the corresponding position. The decreased separation space between fibres in the z direction decreases the loss of the fibre volume fraction caused by the diverging fibres in the y direction (see Fig. 7a.ii). In addition, the fibre orientation in the main members before converging together is slightly oblique towards the boundary in the z direction instead of completely following the topology layout can be observed (see Fig. 7a.iii). This is caused by the feature of the continuity constraint that fibres cannot start or end inside the boundary. To increase the fibre volume fraction, oblique fibre orientations enable extra fibres to join the structure.
This illustrates the constraints for the fibre orientation fulfil the continuity equation and keep the fibre volume fraction as large as possible to minimise the compliance. As a result, a non-uniform fibre stacking in the transverse direction which is not able to be achieved whilst the width and thickness of the deposited fibre bundles are restricted by the nozzle geometry is generated. Despite the continuity constraints encourage the smoothly varying fibre orientation and successfully model the variation of the fibre volume fraction, the manufacturability issue of non-uniform fibre stacking occurs and interpretation of the fibre path is needed.
The second set of the manufacturing constraints considered is the fibre path generation constraints (see Eq. (20)). It is expected to eliminate the non-uniform fibre stacking mentioned before since the aspect ratio and orthogonality constraints are restricting the shape, width and thickness of the deposited fibre bundles. The topology shown in Fig. 6b is thinner in the z direction and wider in the y direction, compared to the continuity case. The fibre orientation comparison in Fig. 7 shows that the fibre orientation at the right of the topology extends straight along x axis instead of diverging along the y axis (see Fig. 7b.i). This shows that the constraints for the aspect ratio and orthogonality discourage non-uniform fibre stacking, however, there may be a loss in terms of mechanical performance due to the fibre orientation no longer being aligned with the principal stress direction. The fibre volume fraction variation stores at the bottom of the diagonal member and the top of the horizontal member where the fibre orientation turns to a closer direction before emerging together at the intersection area (see Fig. 7b.ii). Due to the curvature being larger at the outer side of the turning fibres which leads to an increase in the space between fibres, the fibre volume fraction decreases. This illustrates the decreased fibre volume fraction cannot be compensated by the oblique fibres or converging in the other transverse direction when the manufacturing constraints are considered. From Fig. 7b and Fig. 8a, the streamlines which lie within the isosurfaces of S 1 and S 2 show the cross product of the scalar fields is aligned with the fibre orientation. Overlapping the isosurfaces of S 1 and S 2 gives Fig. 8b where the intersection curves are coincident with the fibre orientation. The cross section shown in Fig. 8b.i where the fibre volume variation is constant shows the isosurfaces form square grids. This illustrates that the orthogonality of the gradients and the aspect ratio constraint are fulfilled. Another cross section shown in Fig. 8b.ii where the fibre volume fraction decreases from the top and bottom to the centre shows the distances between the isosurfaces increase as the fibre volume fraction decreases. This shows how the magnitude of the gradients successfully represent the fibre volume fraction. The cross section at y = 0.5 shown in Fig. 8b.iii also provides an evidence that the distances between isosurfaces are varying with the fibre volume fraction. The intersection curves are then extracted as fibre paths shown as Fig. 8c. The coordinates of the curves give the paths for the print head and the gradients of the scalar fields provide the information about layer height and printing width for slicing software to generate printer toolpaths.
The comparison of compliance shown in Table 4 indicates that the compliance increases as the number of the manufacturing constraints increases. Although the decrease in compliance is 35.69% in the continuity case and 64.65% in the fibre path generation case, compared to the spatially optimised case, the constraints successfully avoid the concerns and problems mentioned in the previous discussion regarding the manufacturability. Despite the compliance of the cases with manufacturing constraints being higher than the spatially optimised case, they are still lower than the unidirectional case. This illustrates that although the variation of the fibre orientation is limited by the manufacturing constraints and the fibre volume fraction is lower at some part of the structure to fulfil them, the design from the present approach has a better performance than the unidirectional case in terms of the compliance when the manufacturability is ensured.

Conclusion
In this paper, a novel design approach to distribute fibres to respect manufacturing constraints has been presented. Using the rule of mixtures to link the fibre volume fraction to the material properties is the first step in connecting the design approach to manufacturing. The desired continuous fibres are ensured, furthermore, when the fibre volume fraction is also coupled to the fibre orientation by the continuity constraints. This brings the optimisation results closer to manufacturable designs. The fibre path generation constraints successfully use the scalar fields to interpret the fibre volume fraction and the fibre orientation which not only reconstructs the model for visualisation but also provides sufficient information for the rest of manufacturing steps.
To demonstrate the potential of fibre-reinforced composites and the effectiveness of the approach, a numerical example with five scenarios is presented. A reduction of 62.35% in the compliance of the spatially optimised case is noted in comparison to the isotropic case. Although the compliance increases 35.69% and 64.65% in the manufacturing constraints cases, compared to the spatially optimised case, the fibre continuity, the variation of the fibre volume fraction and the manufacturability are well controlled in these cases. The comparison to the unidirectional case which is currently manufacturable by commercial desktop printers illustrates that compromising with the manufacturing constraints by decreasing the fibre volume fraction and having the nonoptimal fibre orientation in some part of the structure still results in a better design.
In summary, the present novel approach with manufacturing constraints has shown the ability to generate better designs than existing approaches whilst ensuring the designs are manufacturable.