. Design of 3D and 4D printed continuous fibre composites via an evolutionary algorithm and voxel-based Finite Elements: Application to natural fibre hygromorphs.

This paper describes a novel voxel-based technique to accurately model 3D printed continuous yarns of fibre composites at printed filament scale. To cater for the constraints of continuous filament printing, the design strategy here makes use of a set of variables that fully define the printing path as input in the model. A voxel-based algorithm is developed to obtain the local orientations of the filaments and local quantities of yarns and matrix from a pixelated flat layout of the printed composite. This voxel-based representation is then translated into a Finite Element model to obtain the required structural properties. This approach simplifies and potentially shortens the design of the part, compared to other analytical techniques. The design space of the part is defined by the variables of the gcode used to print the composite. The voxel model is coupled with an evolutionary algorithm to explore the part of the design space related to the use of continuous yarns. The hygro expansion of the material was measured to determine the coefficients of moisture expansion serving as input in the model. Promising results are obtained from a test case related to a continuous flax fibres-reinforced polylactic acid (PLA) printed structure, mimicking the shape of a leaf. The proposed modelling strategy has the potential to characterise geometric, material, mechanical and actuation properties of general 3D printed structures with negligible out-of-plane (i.e., through thickness) printing.


Introduction
Hygromorphs are structures capable of changing shape as a reaction to an external variation of relative humidity (RH) [1].A composite in which each layer deforms in a different manner to the same external stimulus develops a gradient of planar deformations that generate outof-plane bending, hence, an out-of-plane deformation [2].Compared to existing layered multi-material hygromorphs (different material per layer), 3D printing techniques allows varieties of complex architectures, and a broadening of the design space [3].The hygromorph concept proposed by Gladman et al. initiates its deformation at the deposited filament scale [4].Different orientations of the filaments on consecutive layers create hygromorphic deformations due to the hygro-expansion anisotropy of the filament.The flexibility gained by using this material (cellulose-reinforced hydrogel) has however a negative trade-off in terms of poor mechanical properties (longitudinal and transverse stiffness close to 1.27 MPa, while the longitudinal and transverse strengths are 28 kPa and 20 kPa, respectively).4D printed continuous flax fibre reinforced polylactic acid (FFRPLA) composites (printing filament made of a flax yarn covered with PLA, the yarn being made of flax fibres twisted together) are capable of actuation under a Relative Humidity (RH) stimulus due to the flax yarns transverse expansibility [5].These composites also provide load-bearing capabilities in dry (≈11% RH) and wet (≈98% RH) conditions (longitudinal moduli E  ∈ [6.68; 18.45] GPa, transverse ones E  ∈ [0.68; 1.91] GPa, ultimate longitudinal strengths   ∈ [178,201] MPa and transverse ones   ∈ [4.45; 7.97] MPa; values obtained among five RH conditioning: 11%, 33%, 54%, 75% and 98%) [6].
The Timoshenko's model is widely used to predict the changes of curvatures in hygromorph slender beams made of isotropic and anisotropic materials [7][8][9].However, this model is not adequate in the presence of mechanical anisotropy combined with complex stacking https://doi.org/10.1016/j.addma.2022.103144Received 27 April 2022; Received in revised form 19 August 2022; Accepted 8 September 2022 sequences architectures and geometry.A possible way to improve this model is by taking inspiration from stiffness and strength optimisation techniques of 3D printed composites produced with complex printing paths [10].In most cases, the first step to apply successfully these methodologies is to reduce the design space of the geometry of the structure to be modelled.The design space for the geometry could be fairly constrained, for instance, with open hole mechanical test specimens [11][12][13], or extended as it is often the case for truss structures subjected to bending stiffness optimisation [14,15].The second step, specific to anisotropic materials 3D printed filaments, consists in defining the local optimal filament orientations to maximise the required design objective.For instance, Fayazbakhsh et al. have used a sinusoidal function to define the path of their deposited fibres [16].Several research groups have considered analytical approaches to obtain the fibre orientation; for instance, the direction of the principle stress.[11,14], as well as Boddeti et al. who have combined this approach with microstructural homogenisation [15].Other methods, such as bi material discretised optimisation, do not use any analytical equations to distribute the fibre, matrix and fibre orientation through out the structure [17].For certain potential applications (i.e.biomimicking), the main framework of the printing pattern must be defined from the onset.The optimisation in this case serves to tweak the parameters without changing the framework of the printing path.Such strategy lowers drastically the design space, but also makes the optimisation faster due to the reduced number of parameters.Hence, a simplification of the calculations required to acquire the streamlines associated with the material orientation would reduce the computational time and make the implementation easier.The final step is to convert the streamlines obtained for the local direction of the filament to an effective 3D printing path.Kubalak et al. have developed a technique to obtain a printing path from the input of the load path for discontinuous extrusion in 3D and multi-axis printing [18,19].The continuous filament version of this technique has also been implemented [17,20,21].The design space for the printing path of pultruded continuously printed material is, however, restrained, as the filament cannot be cut in the middle of the print prior printing elsewhere.If contact is lost between the filament and the surface, the pultrusion stops.The starting of a print in the middle of a part then generates large defects.For the strength and stiffness optimisations highlighted previously, the objective function is most of the time extracted from the discretisation of the local orientation and the material properties of the printing path [11][12][13]15,16].Finite Element Analysis (FEA) models have also been implemented to describe two-phases [22] and cross-ply continuous fibre composite actuators [23].In the latter case, the hygromorphs have been modelled at the individual layer level.However, in the case of 3D printed composites, the development of a model at the layer scale does not capture the complexity of the fibre (or filament) path involved.Hence, models at the filament scale (the mesostructure) are required to consider the fibre orientation and the amount of material deposited locally in a ply for varying stiffness prints (i.e., structures with varying distance between yarns).Voxel modelling is extensively used in CT-scanning and 3D imaging applications.It is applied, for example, to discriminate the edges of shapes in a given image [24].Voxel modelling has been also implemented to quantify the orientation of fibres in 2D images or 3D volumes obtained from CT-scanning [25][26][27].Voxel-based methodologies have also been used to describe multi-materials structures.Hamel et al. have created a voxel-based FEA model to design hygromorphs with two isotropic homogeneous materials [28].Fayazbakhsh et al. have used a Matlab code to detect defects and the overlaps in automated fibre placed sinusoidal-shaped tape [16].Knowing the size and position of these elements, the properties of the FEA elements are assigned based on the location of the voids (or overlaps), and the elements considered.Those authors propose to assess the potential of processing grey scale images to accomplish similar assignments of the local amount and orientation of the material into the FEA element.
The printing paths of complex geometries are defined by many interdependent parameters.For a given target design, the design space created by these parameters is often too large to be explored manually.Topology optimisation performed via evolutionary algorithms is extensively implemented to optimise the stiffness of complex 3D printed geometries [10].Evolutionary algorithms are able to solve combinatorial optimisation problems in a robust manner [28].For instance, Belhabib et al. have successfully optimised the printing parameters of compression specimens to produce the material with stiffness required for a given weight [29].To the best of the authors' knowledge, such type of evolutionary algorithm has not been so far implemented to design 4D printed continuous fibres reinforced polymers, which are responsive to an environmental stimulus.
The purpose of the study is therefore to develop a novel voxel-based modelling technique that can be used for optimising the filament path of 3D printed composites with continuous filaments, with particular emphasis on humidity-responsive materials.Similarly to what has been presented in open literature, the model has been discretised to run a voxel-based FEA model with varying material properties and the filament orientation of an orthotropic material associated to each pixel.However, the model was adapted to the printing technique used.This technique prevents any cut of the filament during printing.The lack of filament cutting capability limits the design space.To avoid the complex streamline definition for the loading considered, the workflow has been inverted here, compared to what has been described in most of the stiffness and strength optimisations techniques previously highlighted.As such, the input of the model is a set of parameters fully defining a filament printing path.This approach allows to restrain the 3D printing design space to a small part of the continuously printed composites design space.Such restriction is beneficial to the convergence speed of the optimisation.This part of the design space is also defined by authorised values of the set of the parameters fully defining the filament path proposed.This also avoids creating a code that defines the gcode from the streamlines, which is a problematic operation for pultruded materials, as described before.The simplification of the model is made by considering images of the printing path as input in the voxel-model.An evolutionary algorithm is then used to explore the design space.The optimisation technique developed allows restraining the design space of the target component.For instance, the geometries processed during optimisation are all related to the biomimetic leaf design concept used as demonstrator in this work.A continuous flax fibre reinforced polylactic acid hygromorph composite that mimics the structure of a leaf is developed here to assess the potential of this modelling/design strategy workflow.The out-of-plane deformations obtained from this hygromorph are compared to the predicted deformations provided by the model to assess the accuracy of the design strategy.The measurement of the material hygro expansion has been also performed to fully define the material properties required as the input of the FEA model.

Structure of the design methodology
Fig. 1 presents the structure of the code used to design 4D printed hygromorphs.
For the geometry considered in the present study, the smallest Set of Parameters (SoP) fully defining the printing path is depicted in Fig. 2. The design space presents nine dimensions and a potential of 10 9 combinations.
A meta-heuristic method has been chosen and implemented here, because of its simplicity [30].The main issue of such optimisation strategy is the risk of failing to identify the global optimum [30].Due to the large computational resources required for each generation, the optimiser has been set to reach diverse local minima, and not necessarily the global optima.The evolutionary algorithm has been built to make the convergence happening as much as possible on diverse SoPs, in order to limit the risk of producing non printable geometries.To reach this aim, cross-linking and mutating have been made independent from the n 2 best SoP of each step.As described in Fig. 1, n 1 {∖ n 1 ∈ N} SoP are randomly created at the start of the simulation.This ensemble of SoPs is used to make n 1 runs of the voxel model at each generation of the design loop.The output of the voxel run is an ensemble of n 1 fitness values, each corresponding to one input SoP.The fitness of each SoP is compared, ranked and divided in three different subassemblies.The n 2 {∖ n 2 ∈ N} best SoP are selected and goes directly to the next generation, the n 3 {∖ n 3 ∈ N} worst SoPs are completely deleted from the next generation.The rest of the SoPs (the n 1 -n 2 -n 3 SoPs left) go through a function behaving as a genetic algorithm.The SoPs are randomly cross linked (Parameters are swapped between two SoPs) and mutated (Random modification of parameters in randomly selected SoPs).Then, these SoPs are combined with the best sets and completed by randomly created SoPs, to feature the same number of set of parameters than the previous generation (n 1 ).This naturally selecting code iterates until convergence is achieved.The convergence was obtained when the two fitness functions reached 0.0001 mm −1 .The final SoP is checked manually to make sure the corresponding gcode is printable.As a previous study from the research group has shown, the gcode print path sent to the 3D printer is usually not the same as the printed part obtained [5,6].Several local printing patterns are not achievable without deep tailoring of the printing variables, such as printing speed, temperature or layer height.Sharp angles and short straight segments (2-3 mm) are complex to print, due to the reduced pulling force from the support surface to the pultruded filament.Hence, each print path has to be adapted to obtain the exact print path obtained from optimisation (i.e. for short segments they must be made longer).This adaptation of the printing path depends on the printing technique, the material considered and the experience of the operator.The geometry is then exported as a gcode to be printed.The algorithm is run on the High Performance Computer Blue Crystal3 (University of Bristol), under a Linux operating system, 57 GB of RAM and 16 processors on one node.Fig. 1 shows that the geometry to be designed must be defined first, before launching the design algorithm.
A fitness function is used to define the best and the worst SoPs in one run.This objective optimisation function () uses the difference of curvatures between dry and wet states at specific locations of the specimen undergoing the hygro-morphing.The function can include as many curvatures as needed: n {∖ n ∈ N}, and it is in Eq. ( 1), in which k  is the curvature obtained from the FEA model and kref  is the curvature target.Finally, a  {∖ a  ∈ R} are coefficients defining the priority convergence between the curvatures: Additional constrains can be applied for the selection of kref  to avoid degenerate solutions, in which curvatures set as targets are outside the achievable range, or not compatible between them.The range of curvature achievable and the compatibility with local curvature depend on the material implemented, the conditioning RH, the geometry and filament path considered (In the present paper, ''conditioning'' refers to immersion in a environment with a controlled relative humidity (RH).Experimental or preliminary evaluations via model or experiment are required to define these limits and the above compatibility.The evolutionary algorithm also serves as a way of understanding if the local curvatures used as target can fit together.The Supplementary material provides an example of degenerated outcome for the optimisation process.
The optimisation algorithm requires solving the following problem: minimise  subject to Equation (1) The term  indicates the interpolation function that defines the curvature from the displacement of the nodes   .K and F ℎ are the global stiffness matrix and generalised force vector associated to the change of relative humidity in the environment around the hygromorph, respectively.The terms  and   indicate the orientation and the amount of the material (fibre volume fraction) in a given element of the FEA model, respectively.The discretisation of the problem implies to integrate the function F that defines the orientation and the fibre volume fraction in a given set of SoP  over the nine-dimensional design space .The SoPs are distributed between the upper (  ) and lower (  ) bounds, the latter defined after assessing the real capabilities of the printing.
For the geometry chosen as an example in this study, the specifications include the size of the leaf to be circa [120; 130] mm (l 2 ) long and [25; 30] mm (l 1 ) large.The target curvatures are 0.0094 mm −1 along the length direction (k 2 ), and 0.007 mm −1 on the transverse direction (k 1 ).These specifications are shown in Fig. 3.The optimisation algorithm is terminated as soon as the two fitness measured reach below 0.0001 mm −1 .

Voxel-based modelling
approach [24].The eigenvector corresponds to the smallest eigenvalue of the Hessian matrix of the pixel analysed (Example displayed in Fig. 4(b)).The Hessian matrix (H) is determined via Eq.( 3).The fibre volume fraction associated to each pixel is estimated from Eq. ( 4) (Example displayed in Fig. 4(c)).
where f is the grey scale image of the filament printing path.
The term BW indicates the shade of grey on a black and white scale.
The component   corresponds to black level on the black and white scale.The resolution of the image is tailored to obtain the adequate number of pixels for the width of the yarn.The following steps have been implemented to obtain the adequate dot per inch (dpi).• When the .inpfile is created, the dimension of the model must be checked.The amount of pixel must be adjusted to the computing capability available.
Each pixel corresponds to a FEA element with specific orthotropic material properties and orientation.The S4RT shell element in ABAQUS (Dassault Systems) is used in the present model.The humidity-induced deformation is simulated via thermal boundary conditions loading made available in Abaqus.The model built for the leaf test case conducted between 11% and 98% RH.Such levels of humidity lead to a 5.77 * 10 −2 % moisture content variation (loading used in the FEA) inside the flax fibres reinforced PLA (FFRPLA) composite (See the hygro expansion results in Section 5.4).The out-of-plane displacement and the displacement along the length direction (z and  in Fig. 3, respectively.) of the entire sample of the specimen are fixed at its triangle shaped tip (point (b  ,b  ) in Fig. 2).These displacement constraints have been observed to provide a small impact on the curvature obtained from the specimen.The properties implemented in the model are presented in Table 1.To choose the properties to be implemented, priority has been given to the characteristics of the FFRPLA composite measured via tensile and hygro expansion tests made in the present study, and in a previous work performed by the same authors [6].In this study, the properties required have been complemented with engineering constants  23 ,  23 from open literature.For  23 , the value related to flax fibres was considered.The Chamis micromechanical model has been used to obtain the equivalent composite transverse shear modulus  23 [31].In this equation, the shear modulus of the PLA is considered equal to G   = 552.6MPa, and the fibre volume fraction is 0.3 [6].The shear stiffness (G   ) of the PLA matrix has been determined via the relation between tensile and shear modulus for isotropic materials    = 0.36 [32], with E   = 1503 MPa [33].
All material properties (except for the Poisson's ratio) are multiplied by the fibre volume fraction (FVF) of the corresponding pixel.All pixels with FVF=0 have their material properties multiplied by 10 −6 , to avoid numerical instabilities compared to the properties of pixels with FVF≠0.

Hygro expansion and coefficient of moisture expansion
The hygromorph 3D/4D printed composites can expand along the three dimensions due to moisture effects.The specimens used for the measurement of the hygroscopic expansion are 20 mm long and wide, and 1.2 mm thick.These specimens consist of five layers of 0 • orientated flax yarn reinforced PLA.The deformation is measured via a Mitutoyo micrometre IP65 with a 0.001 mm resolution.The Fisher Scientific PAS214C scale used to determine the mass has a 0.0001 g precision.The drying of the samples is performed at 40 • C under vacuum for at least 72 h.After drying, the specimens are placed in hermetically closed boxes containing de-ionised water saturated with different salts to obtain specific Relative Humidity (RH) values.The duration of the environmental conditioning varies according to the different salts used.The specimens have been tested when their weight has reached stability after conditioning.The RH values used in this work, the conditioning time, together with the types of salts adopted for the conditioning of the samples are presented in Table 2.
The moisture content (MC) is calculated using the mass of the material in dry state   and the mass at the considered RH value   : The coefficient of moisture expansion (CME -) is measured along each dimension via the following Equation: In Eq. ( 6), l  and l  are the dimensions of the dried and conditioned specimens, respectively.

Printing
The printing of continuous FFRPLA composites has been described in a previous study [35].The printing path is defined using a Matlab code and it is used to evaluate the precision of the implemented FEA model.Fig. 5 shows the printed geometry, as well as the parameters implemented for the printing of the test case considered throughout this study.The optimisation algorithm shown in Fig. 1 specifies that the optimised printing path must be first checked against the real printing capabilities of the additive manufacturing machine.In this case, the printing path related to the specific printed geometry must be modified to be as close as possible to the optimised configuration.For the geometry framework presented in Fig. 2, the value of   has first been checked: if too large, the filament would not be able to stick to its support layer.All lines before and after a back and forth corner have been increased by 1.2 mm in length, to allow the filament to stick to its support layer.This phenomena was explained in a previous study of the research group [6].

Validation of modelling and hygromorph simulation
In the first conditioning step, the printed specimens are subjected together to 98% RH for 12 h.The specimens are then oven dried one by one at 38 • C for 6 h.While one specimen is oven dried the others are placed in the box at 11% RH.The samples are then transferred for at least five and a half days in a sealed box containing water saturated with potassium hydroxide, lowering the RH of the box to 11%.The specimens are then transferred to a desiccator with a 98% RH inside.The 98% humidity is obtained after a preliminary study in which crystallisers filled with distilled water and others filled with potassium sulphate saturated distilled water were positioned in the desiccator used in the experiment.A Testo 174 h sensor has been used to monitor the correct range of humidity.The specimens deformations are measured via a 3 Dimensional Digital Image Correlation (3D-DIC), and using a Davis Lavision 5 Mega Pixel camera set up with 29 pixels subset size and 9 pixels step size.Sum of differential correlation mode was implemented for the postprocessing.The specimens have been manually speckle patterned to control the size of the points.The samples were glued to a transparent and translucent string to be held vertically.The vertical position of the specimens limits the friction during deformation, but also avoids a differential moisture content through the thickness of the composite itself (Fig. 6).As the desiccator is closed, the humidity raises from the bottom of the box; if the specimen was horizontal, the humidity would reach first the bottom surface and pass through the thickness later.Therefore, the resulting transition state of the deformation would not be representative of the uniform relative humidity intake simulated in the FEA model.
The DIC records one image every 30 s for the two first hours of deformation, followed then by one image every 5 min for a week.The post-processing of the deformation measurement is presented in Fig. 7.The DIC subsets have been plotted in 3D, then an interpolation function consisting in a fourth order polynomial along two spatial variables ((x, (y) ⧵ P(x,y)=z) has been considered for the fitting of the geometry.Fourteen points along the axis of the curvature are tracked during the deformation and the duration of the test.The values used for the design comparison are determined with Eq. (7).Examples of curvatures extracted from this type of calculation are given in Fig. 7.
=  − −  − (7) In Eq. ( 7), i indicates the index of the curvature defined in Eq. ( 1).The terms a and b are the indexes for the initial and final deformation stage, respectively.

Validity of the FEA model
The validity of the FEA model is assessed via comparison with the printed structure.Fig. 8 shows the evolution of the curvature during the deformation measurements performed via DIC.The two curvatures used in the optimisation function (see Fig. 3) are presented in these graphs.Both curvatures converge over one week or more.Specimens for a given curvature converge at different pace.An end criterion for the test had to be defined for the comparison of the same quantities obtained by FEA.The test was considered terminated as soon as the curvature varied less than 3% in 48 h.The experiments ran up to 21 days, depending on the specimen.The curvature obtained at this point (3%, 48 h) was also used for comparison with the FEA model (Fig. 8).The curve for the evolution of k 2 are interpolated with double power functions (ki 2 (t)=a1*t 1 + a2*t 2 with {(a1,a2,b1,b2)∈ R 4 }).All these interpolations presenting an average correlation coefficient of 0.98 ± 0.02, the specimen could be concluded to have similar deformation capability.On the other hand k 1 does not present such clear similitude between the different specimens.Interpolation with an exponential function (ki 1 (t) = a1*e 1 *  + a2*e 2 *  with {(a1,a2,b1,b2)∈ R 4 }) was implemented after the maximum of the curve was reached.However, the diversity in curvature evolution led to an average correlation coefficient of 0.86 ± 0.10.The standard deviation to average ratio of the curvature measured at comparison point presented in Fig. 8 is 60% for k 1 and 44% for k 2 .The large discrepancy is also explained by printing difficulty.The leaf-inspired hygromorphs printed in the present study are, to the author's knowledge, the first structure with such complex filament path to be printed out of 3D printed continuous flax yarn reinforced PLA.There is still significant room for improvement regarding the precision of the filament deposition.Fig. 5 shows several features in the specimens arising from printing issues, which highlight the difference between the geometry modelled and the ones printed.The curved print of the upper tip of the specimen does not entirely sits on the straight filament layer, and this is an issue mainly due to the difficulty to print short segments with the current rig, as the filament does not perfectly sticks to its support layer.However, these printing issues are very complex to account for and varies between specimens printed with one gcode.The variability of the flax-PLA filament (amount of resin, natural origin of the fibres, twisting of the yarn) used has an impact on the variability of the print obtained [35].As previously discussed, the variability in the print and the difficulty to print certain patterns is the first reason for the discrepancies between curvature values measured within the different specimens.The variability of the curvature measured between specimens also comes from the measurement technique used in this work.The rig does not allow to measure accurately deformations over small distances.The ratio between the standard deviation and the average values is almost twice as high for  1 than for k 2 .The order of magnitude for the two curvatures is identical.However, k 1 was measured over 28 mm, whereas k 2 was measured along a 124 mm straight line.The longer the distance used to measure the curvature, the greater the amplitude of the in-depth position; the interpolation of the curvature will be therefore less sensitive to the variability of the position.Zhang et al. modelled via FEA 3D printed continuous carbon filament to observe the influence of the printing path on the stresses inside the filament as well as on the geometric features of the print [36].Those authors highlighted that the reduction of the printing radius in corners creates many defects, such as twist of the filament, error path, fibre breakage, folding and misalignment of the filament itself.Zhang et al. have proposed a FEA model at the filament scale as an efficient approach to predict those defects.However, the scale they have implemented cannot be considered to model large structures like the one implemented in the present study, also because of the limitations of the measurement technique.The initial stage corresponds to the first image recorded, as the specimen is positioned in the measurement rig presented in Fig. 6 (11% RH conditioned specimen).The final stage corresponds to the moment the test was considered as terminated (i.e., the specimen conditioned at 98% RH).The termination criterion is explained in the Discussion section.The position coordinates (x, y, z) refer to the coordinate system shown in Fig. 6. between filament the larger the deformation capability of the yarn (confirmed for flax fibre reinforced polymer [37]).The greater the compaction the smaller the actuation.
Fig. 2 shows that the compaction (distance between yarns) of neighbouring filament varies significantly from the position between layer but also through each layer.Ideally speaking, each voxel would have had its CME modified depending on the impact of its neighbouring yarns.However, such complex characterisation is not achievable for now.Consequently, the CME considered is chosen to avoid too constrained or loose hygroscopic deformations.To this extend, the through-the-thickness CME presented in Fig. 12 was considered for the model.On one hand, the in-plane coefficient of moisture expansion (  = 0.123) presented in Fig. 12 is adequate to model the triangular tip of the geometry printed in Fig. 2, as this part of the printed structure present filaments compacted against each other with a 0.8 mm distance between successive filaments.However, the in-plane CME along the width is not capable of reproducing the behaviour of the rest of the structure, due to the lower compaction of the filaments.That CME value artificially stiffens the modelled structure and leads to model deformations (k 1 = 0.0031 mm −1 ; k 2 = 0.0025 mm −1 ) that are significantly lower than those measured.On the other hand, the values related to flax fibres present in open literature (  = −0.00629[38];   = 1.060 average from [38][39][40]) do not account at all for the constraining of the neighbouring filaments.Those values artificially soften the modelled structure, leading to larger curvatures (k 1 = 0.026 mm −1 ; k 2 = 0.013 mm −1 ) than those measured in the experiment.Part of the discrepancy between models and tests is also coming from the large difference observed between the different tested specimens, as explained in the previous paragraph.The in-plane motion of yarns and fibres reinforcing hygromorphs has been highlighted by FEA model in open literature [41].Similar deformations are happening in the structure presented in the paper, introducing differences between the model and the test case adopted here.Finally, the printed geometry also does not exactly match the theoretical target printing path.Unlike the printing pattern presented in Fig. 2, the leaf-like structure features sharp corners that are extremely problematic to obtain when the direction of the filament changes abruptly (see Fig. 5).The difference between the printing path and the printed geometry at the corner is small, however the presence of a large number of corners affects the overall deformation of the printed specimen.

The optimisation algorithm
Fig. 10 shows an example of the evolution of the best fitness during the design run of a leaf structure.This figure also presents the evolution of the standard deviation of the fitness function for the n 2 (see Fig. 1) best SoP.The tracking of this parameter allows to ascertain if the different optimised SoPs are usable to design the required geometry.The identification of several possible geometries adapted to production is made, in case the best one is practically identified as not printable.The experienced 3D printing user can limit the number of non-printable geometries by reducing the interval of the selection of the values that define the specific geometry.For the algorithm charted in Fig. 10, the convergence coefficients are 0.1 and 0.9 (respectively, a 1 and a 2 from Eq. ( 1)) for k 1 and k 2 , respectively.The fitness of k 2 converges fast towards a null value; it is therefore possible to prioritise certain curvatures within the design of a particular class of hygromorphs.The optimisation technique was demonstrated as capable of converging towards new SoPs for curvatures known to be achievable prior to optimisation.
The processing time and CPU required for the convergence are mostly dependent on the size of the voxel model.That is why the choice of the pixel size is fundamental, and must be limited as much as possible when building the code.Among the trials carried out for the design optimisation of the leaf-inspired hygromorph a large range of convergence timing were obtained.The example presented in Fig. 10 converged in three days of computations on Blue Crystal 3. Others did not converge in regard to the criteria set up despite a 3 weeks run.However, studying the impact of the workflow parameters (n 1 , n 2 n 3 ...), the target curvature and the model input properties is not part of the discussion of this study.Fig. 11 shows the distribution of the design parameters for the optimised geometry at the end of the run of the optimisation algorithm.
The large scatter of the parameters considered in this work show that the algorithm does not converge towards one specific geometry.The large number of SoPs deleted at each generation (n 3 :20 out of n 1 :42) allows the algorithm to explore several locations in the design space.The broad sets of SoPs evaluated by the algorithms also allow the identification of local minima from diverse regions of the design space.An issue intrinsic to this methodology is the characterisation of these local minima.For instance, the fitness of the n 2 parameter could be locally improved because that parameter is not involved in the cross-linking and mutation process.In the case of use of a less computationally expensive algorithm, or of a less constraining printing method, it is advisable to include the n 2 best SoPs in the cross-linking to explore to a greater extend each local minima, and potentially find the global minima.Fig. 11 shows that a large part of the available interval is unused.If another run for the same target geometry is to be conducted, these unused part of the parameters interval could be deleted.This would enable faster convergence and greater exploration of a more restrained part of the design space.

Discussion around the overall design strategy
One can identify few advantages and weaknesses of the design strategy and the modelling technique implemented in this work to model structures and hygromorphs produced with continuous 3D printed filaments.detect the contrast of the space between yarns.This is essential to position efficiently the yarn in the FEA model.Fully-fused actuating materials would be required to adapt the Hessian matrix calculation.

Hygro expansion properties
Fig. 12(a) shows the evolution of the Moisture Content (MC) inside the flax fibre reinforced PLA with the change of the Relative Humidity (RH).A mathematical model for the variation of the MC in flax fibres versus the surrounding relative humidity has been previously described by Hill et al. [42].Fig. 12(b) shows the impact of the MC on the Hygroscopic Strain of the specimens.The variation of the length is null at any MC level.The sigmoid function has been previously proposed to interpolate the hygroscopic strain [23].As described in Fig. 12, the sigmoid is an efficient function to interpolate width and thickness.The CME for the specimens between dried and 98% RH are found to be 0.394 ± 0.055 and 0.123 ± 0.010 along the thickness and width, respectively.The CME expansion along the length direction is two orders of magnitude smaller than along the width and transverse directions (−0.0071 ± 0.005).Those coefficients are to be considered as linear approximations of the moisture expansion effects, because they are extrapolated from classical moisture expansion relations (see Eq. ( 6)).

Conclusion
of the design space allowed to shorten convergence.Nonetheless, the modelling and optimisation technique described in this work could also be extended to 3D/4D printed composite structures using different classes of long filament reinforcements.Fossil and non-fossil long fibres could be considered and implemented in the design of adaptive smart 3D printed composites.Homogenisation of local repetitive patterns is one of the options considered to improve the computational efficiency of the model.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Workflow of the design algorithm (SoP: Set of Parameters).Possible values of each parameter in the SoP are defined in continuous intervals.These intervals are discretised to limit the number of possible geometries.The optimisation function is defined in coming paragraphs.Additional constraints to be implemented for the geometry can be also defined (i.e.general dimensions of the print).n 1 : Number of SoP in the population; n 2 : Number of SoP selected at each generation to go directly to the next generation; n 3 : Number of SoP deleted at each generation.

Fig. 2 .Fig. 3 .
Fig. 2. Mesostructure and definition of the parameters defining the printing path for the leaf-inspired hygromorph used as model validation example: [(y); ; n  ; g .10 ; g .11 ; g .20 ; g .21 ; g .1 ; g .2 ] (a): General geometry of the specimen (Not to scale).The coordinate systems refers to the one presented Fig. 3. (b): Image of a leaf: The blue arrows points to the secondary vein crossing the width of the leaf.(c): Printing path for the first layer; polylines inspired by the secondary vein structure of a leaf.(d): Printing path for the second layer; quadratic functions defined to produce as much local cross-ply composite as possible [0 • ,90] to obtain out-of-plane deformation.Definition of parameters : function to define the angle of the filament path depending on its position along the -axis.: angle defining the shape of the triangular tip of the print.n  : number of paths printed in the second layer intersecting the vertical axis of symmetry (10 in (d)).g: distance between successive filaments; g  and g  refer to first (transverse) and second (longitudinal) layer respectively; The following indices (i;j) in g . and g . i corresponds in the section of the print (i = 1: triangular section; i = 2: curved section), j is used to distinguish between the different interfilament distances used over this section of the print.(e): Example of deformation of the leaf structure observed with the digital image correlation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. General structure of the voxel code.This flow chart is inserted in the place of the Voxel model box in Fig. 1. (a): Images of the printed path.Greater details on how to obtain this image is presented in supplementary material.(b): Each pixel has the orientation of the curve calculated and displayed.(c): The amount of composite for each pixel is displayed in these figures.

Fig. 5 .
Fig. 5. Printed geometry for the SoP presented on the right.The definition of each parameter in the SoP has been presented in Fig. 2 This geometry serves for the comparison with the voxel model.(a) Layer 1: bottom up view, (b) Layer 2: top down view.

Fig. 6 .
Fig. 6.Schematic of the set up used to record the deformation of the leaf inspired actuator.Two different views are presented: (a) Front view, (b) Top-down view.

Fig. 7 .
Fig. 7. Evolution of the leaf-inspired shape before and after conditioning.Measurements made via Digital Image Correlation and presented in different views: (a) lateral, (b) top down and (c) isometric view.The terms k 1 and k 2 refer to the curvature presented in Fig. 3.The subscripts a and b refer to the initial and final deformation stage, respectively.The initial stage corresponds to the first image recorded, as the specimen is positioned in the measurement rig presented in Fig.6(11% RH conditioned specimen).The final stage corresponds to the moment the test was considered as terminated (i.e., the specimen conditioned at 98% RH).The termination criterion is explained in the Discussion section.The position coordinates (x, y, z) refer to the coordinate system shown in Fig.6.

Fig. 8 .
Fig. 8.Comparison between the evolution of the curvatures predicted by FEA and measured via DIC in specimens conditioned at 98% RH, and related to the geometry presented in Fig. 5: (a) for k 1 and (b) for k 2 .Each blue shade curve present one of the four specimens tested.The two horizontal dashed lines present the deformation obtained via FEA model.The purple error bars introduce the average and standard deviation of the measured curvature at the point of comparison between the FEA model and the test as defines previously; (3%, 48 h).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .Fig. 11 .
Fig. 10.Example of fitness convergence obtained for the leaf design.The black markers are indicative of the best stiffness at each generation.The blue and turquoise markers present the fitness term of Eq. (1) related to k 1 and k 2 , respectively, for the best general fitness.The error bars show the standard deviations for n 2 values of ,  1 and  2 (Fig. 1).The parameters used for the run are: n 1 = 42, n 2 = 4 and n 3 = 20.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2
Relative humidity values used for the conditioning in climatic chambers and the corresponding salts used.An estimation of the duration of conditioning to obtain stable shape is presented in the last row.