Optimal design and manufacture of variable stiffness laminated continuous fiber reinforced composites

Advanced manufacturing methods like multi-material additive manufacturing are enabling realization of multiscale materials with intricate spatially varying microstructures and thus, material properties. This blurs the boundary between material and structure, paving the way to lighter, stiffer, and stronger structures. Taking advantage of these tunable multiscale materials warrants development of novel design methods that effectively marry the concepts of material and structure. We propose such a design to manufacture workflow and demonstrate it with laminated continuous fiber-reinforced composites that possess variable stiffness enabled by spatially varying microstructure. This contrasts with traditional fiber-reinforced composites which typically have a fixed, homogenous microstructure and thus constant stiffness. The proposed workflow includes three steps: (1) Design automation—efficient synthesis of an optimized multiscale design with microstructure homogenization enabling efficiency, (2) Material compilation—interpretation of the homogenized design lacking specificity in microstructural detail to a manufacturable structure, and (3) Digital manufacturing—automated manufacture of the compiled structure. We adapted multiscale topology optimization, a mesh parametrization-based algorithm and voxel-based multimaterial jetting for these three steps, respectively. We demonstrated that our workflow can be applied to arbitrary 2D or 3D surfaces. We validated the complete workflow with experiments on two simple planar structures; the results agree reasonably well with simulations.


Anisotropic Filtering
Figure S1 Illustration of anisotropic filtering for the hyperbolic paraboloid shaped laminate. The blue oblate spheroid represents the filtered region. The solid arrows indicate the local coordinate system used to orient the filtered region and calculate the anisotropic filter weights while the dashed arrows indicate the global coordinate system.
Filtering of the optimization variables is known to help regularize the optimization problem [3][4][5] . The most commonly used filtering approaches include linear 4 and projection filters 5 . Linear filter, in the finite element setting, averages variables defined on FE mesh nodes (optimization variables, , in our case) within a specified neighborhood, e.g., a sphere or ellipsoid. This helps smooth out undesired large gradients in the design parameters (density, fiber volume fraction and orientation). The projection filter is applied to mitigate the blurring introduced by the linear filter and to obtain a distinct 0-1 design, i.e., no intermediate values making it easier to interpret the resulting design. Projection filter achieves this by mapping the input, typically defined over [0,1], to either 0 or 1 with the help of mathematical functions like smooth-Heaviside step or sigmoid. The linear filter is combined with a projection filter to achieve a minimum length scale when applied to density, or volume fraction, . 4,5 In this paper, each design parameter was defined as an elemental property, = ( ) ( could be , or ), using a filter, operating on the nodal optimization variables, ( ⊆ ). This filter operator was set to be either a linear filter, ( ) or a combination of linear and projection filters, ( ( )). Thus: Here, and are the bounds set on the design parameters (see Table S1). The linear filter operation is defined as a weighted average of : Here, is a matrix containing weights assigned to neighboring . The parameters and are the radii of an oblate spheroid (ellipsoid with two of three radii equal) that define the neighborhood over which is averaged. We set = for and , making it a spherical neighborhood, but use different values for to decouple ( >~ layer thickness) the influence of optimization variables across different layers of the laminate. To do this, we transformed the ellipsoidal neighborhood at the point of interest such that the axis with smallest radius, gets aligned with the local normal of the surface as shown in Fig. S1. The elements of the matrix were defined as: Here, are the weights assigned to the pair of i th element and j th node. These depend on a distance measure denoted by and defined as: Here { , , } are the coordinates of the centroid of element i while { , , } are the coordinates of node j. It is to be noted that the linear filter averages = − /2 and = /2 to 0 even as = ± /2 result in the same material stiffness tensor.
The projection filter is used to demarcate the different material phases used in the optimization problem. Typically, with isotropic materials, it is applied on linearly filtered . In our case, we applied it on either or as needed. The projection operator is defined as: Here, and are parameters that control the sharpness and threshold of the filter, respectively. Based on our bounds for , varies between and 1, i.e., 0 ≤ ≤ ≤ 1 where = 0.001 for and 0 for . Values of below are mapped to 0 and above to 1. To avoid numerical issues related to the large gradients associated with this filter, we increased from 0 to 50 and decrease from 0.5 to 0.1 every 100 iterations in four steps. As stated in the main text, the objective in all the topology optimization problems is to maximize stiffness. This was achieved by minimizing the resulting strain energy in the structure due to the applied load and boundary conditions. A constraint on either overall material usage or fibrous material usage was imposed. Tables S1 tabulates the specific constraint used for each problem along with design parameter bounds and anisotropic filter parameters. In Table S2, we list the symmetry conditions imposed on each design parameter.

Simulations: Setup
These symmetry conditions were imposed explicitly for the cylinder and hyperbolic paraboloid problems where a full model was used for FE simulations. For the other problems, we used only a quarter of the design domain for FE simulations while assuming a symmetric structural response. Thus, most of the symmetry conditions were imposed implicitly. The symmetry about the diagonal plane for the plate and spherical shell problems were imposed explicitly. For all the laminated structures considered, we used hexahedral solid shell elements based on Schwarze and Reese's work 6 .

Simulations: Thickness Optimization
The anisotropic filtering scheme enables us to optimize thickness of the laminate along with the microstructure. To demonstrate this, we set up an identical problem to the laminated plate optimization problem (see Fig. 3, main text) with only one difference. We modified the linear filter applied to densities to an anisotropic one with = 0.4 mm and = 0.1 mm. The reduced filter radii have the effect of reducing the smoothening effect of the linear filter and enables appearance of thinner features locally, more so in the thickness direction. This manifests as an ability to optimize thickness as can be seen in Fig. S2. It is to be noted that the resulting thickness variation is continuous and hence results in thin sections where fibers cannot be placed with the fabrication technology we used. Hence, the design with isotropic density filtering was preferred in the current study.  Fig. S4, while Figs. S7-S9 show the same for nonplanar problems. For the problems where the projection filter was employed, the effect of the change in the filter parameters, and , can be clearly seen as spikes in the constraint graph (or dips in the objective graph) at 100 th , 200 th , 300 th and 400 th steps (see plots on the right in Fig. S3 and Figs. S7-S9). In all the optimization problems, we normalized the objective with the strain energy of the initial design and chose an initial design with a uniform or that satisfied the constraint. The initial fiber orientations were, however, chosen through trial and error to avoid local minima. We considered the optimization problem to have converged when the constraint was satisfied and the change in the objective value between two consecutive design iterations stays below 1e-4 for ten consecutive iterations. In all the cases, we ran the optimization procedure past the point of convergence until the algorithm reached a set maximum number of iterations (400 for problems without projection and 800 with projection). Interactive three-dimensional visualizations of the designs are available on the author's webpage: https://m3lab.github.io/cont_fibers.html.

Material Compilation
We present here the key mathematical details of Knöppel et al.'s 7 stripe patterns algorithm. For a complete and detailed treatment, the readers are referred to Knöppel et al.'s 7 paper. To reiterate from the main text, the goal of the stripe patterns algorithm is to obtain a global, periodic scalar field or parameterization, aligned with an input unit vector field . Mathematically, this is equivalent to finding a scalar field, such that its gradient, ∇ = ⊥ , with being the prescribed frequency of the stripe periodicity and ⋅ ⊥ = 0. Assuming = arg( ) ( is a complex valued function with a unit L 2 -norm, i.e., ‖ ‖ =1), this condition is transformed to: Use of the arg function ensures the periodicity of . In the discrete setting, Eq. S6 can be easily integrated along an edge of the mesh elements ( = edge vector connecting vertices and ): The LHS of Eq. S7 can be simplified to ∫ ∇ln( ) ⋅ = ln ( ) by way of the fundamental theorem of calculus for line integrals. Denoting, = ⊥ and assuming varies linearly along the edge simplifies RHS in Eq. S7 to 1 2 ( ⋅ + ⋅ ) ≔ . Thus: Eq. S8 suggests that the argument of needs to change by going from vertex to for the parametrization to align with the input vector field. Since this condition cannot be satisfied exactly for all the edges, the L 2 error, ( ) is minimized. This is defined as: Here, are edge weights that account for the geometry of the mesh elements. For triangular meshes, = 1 2 (cot( ) + cot( )) with and being the angles opposite to the edge within the two triangles that share the edge. For boundary edges, the unknown angle is set to /2. Necessary adjustments are made for edges and triangles with non-orientable and locally orientable (see Section 3.3 of Knöppel et al. 7 ).
The energy ( ) can be represented in matrix form as, ( ) = x x where x is a column matrix populated with the real and complex parts of . Likewise, the unit norm condition ‖ ‖ =1 can be discretized as x x.
Matrix is symmetric positive semi-definite and has structure similar to discrete mesh Laplacian matrix, while is a diagonal matrix whose diagonal elements depend on the areas of the triangle (see Section 4.1 of Knöppel et al. 7 ). Thus, in the discrete setting, the problem distills to: . . = 1 (S10) This is equivalent to solving for the smallest eigenvalue, and the associated eigenvector of the generalized eigensystem, x = λBx.
The parametrization, is obtained from the solution of Eq. S10 by first extracting .Using arg directly to obtain at each vertex will lead to severe aliasing on coarse meshes. Thus for each triangle , is calculated at one arbitrarily chosen corner as = arg and values at the other two corners are locally adjusted using and optimal (see Section 4.3 of Knoppel et al). Linear interpolation of over each triangle and use of a periodic color function i.e., ascribing the same color to values separated by 2 ( is an integer), then results in density plots as shown in Fig. S10b, S10d and S11b. Iso-2 contours then give the desired fiber layout as shown in Fig. S11c. Figure S10 demonstrates the stripe patterns algorithm with two simple examples. In the first case, the input unit vector field is 1 = ( , − ) where 2 = ( 2 + 2 ) with and being the x-and y-coordinates respectively. The perpendicular vector field, ⊥ 1 = ( , ) has zero curl and thus results in stripe patterns with no singularities except near the origin where the input vector is not well defined (see Fig. S10b). In the second case, the input vector field is 2 = ⊥ 1 . The curl of the perpendicular vector field is non-zero everywhere which leads to several singularities as can be seen in Fig. S10d.  Figure S11 illustrates the material compilation workflow we adapted in this paper. The workflow begins with the stripe patterns algorithm which takes as input the optimized design obtained from the design automation step as an input. To be more specific, the input is a surface mesh with spatially varying stripe frequency, and input vector field, defined at each nodal point. The surface mesh is the mid-plane of the shell elements of each layer used in the FE model (e.g., 8 layers for the laminated plate example, 2 layers for the cylindrical shell). The optimized fiber orientations, , on the laminate surface were used to specify the vector field while stripe frequency is defined in terms of the optimized fiber volume fractions, . To establish the relationship between and , we assumed that the fibers fabricated via material jetting are cylindrical with radius and are enclosed in imaginary concentric cylindrical tubes that are hexagonally close packed when viewed along the cross-section perpendicular to the local fiber direction. Thus, the diameter of these enclosing cylinders, , gives us the spacing between each fiber and its six equidistant neighbors (assuming uniform packing). The ratio of the cross-sectional areas of the enclosing cylinders and the fiber should then equal their packing densities. The packing density, , for the enclosing cylinders is simply given by = /√12and the packing density for the fiber is . This then gives us a relation between the fiber spacing, and volume fraction: . Stripe frequency is then simply = 2 / . These relationships were found to extend well to non-uniform .
It is to be noted that the fiber orientations can be flipped by 180° on the laminate surface without affecting the FE simulations. We essentially chose the direction of the vector field, , to be arbitrary and found it to not affect the stripe patterns and thus the fiber layout significantly. The output of the stripe patterns algorithm is shown in Fig. S11b, where the colors represent the scalar field, . It can be clearly seen that the stripe frequency correlates well with the input volume fractions. In Fig. S11b, we also illustrate the ability to use dissimilar meshes for FE simulations and material compilation. This lets us use fine and coarse meshes for FE and material compilation respectively, thus making the material compilation step fast and computationally efficient. We used extrapolation functions to define and in the non-design domain regions and trimmed them away after the fiber layout was generated as can be seen in Fig. S11c.

d) the 3D voxel model containing the voxelized fibers (magenta) and matrix (light gray) to be used for voxel-based 3D printing
The fiber layout, as shown in Fig. S11c with a single layer of fibers, was obtained with the help of VTK Python wrappings. VTK (visualization toolkit) is a scientific visualization software with several algorithms implemented as filters that act on input datasets to output a modified dataset. These datasets could include a set of points, lines, polygons and polyhedra. In our case, the input dataset is a surface mesh with the calculated scalar field from the stripe patterns algorithm stored as a texture coordinate. The fiber layout was extracted as iso-contours of the texture coordinate field with the help of vtkContourFilter. There is no unique fiber layout, we just choose one among several different sets of contours with the same and distributions that could be obtained from the texture coordinates. For example, the fiber layout in Fig. S11c and Fig. 2d of the main text have the same fiber orientation and volume fraction distribution. The layout in Fig. S11c was obtained by simply offsetting the iso-contour values by half of the fiber spacing value. Other VTK filters were employed to stack the layers of fibers, clean up and save the resulting fiber layout for voxelization. To enable manufacturing with voxel-based 3D printing techniques, a voxel model of the structure to be printed is created with a voxelization routine. Figure S11d shows the top view of the voxel model for a single layer of the composite structure.
A summary of the key steps involved in our material compilation algorithm follows. For each layer of the optimized laminated structure: 1. Prepare a triangulated surface mesh based on the optimized design.
2. Store at each vertex a vector, = / ( is a unit vector along the optimized fiber orientation). Extrapolate values where needed. 3. Run stripe patterns algorithm with the prepared mesh and desired fiber spacing as inputs a. Extract edge data: b. Construct energy and mass matrices: and c. Solve for the smallest eigenvalue of the eigensystem: x = x d. Extract the parametrization, from x 4. Extract iso-2 (or iso-2 + , is any arbitrary angle) contours

Algorithm for Connected Fiber Layout
We enumerate here the steps involved in the algorithm for obtaining a connected fiber layout. The input to this algorithm is the disconnected fiber layout comprised of polylines as obtained from the material compilation step.
1. Identify boundary vertices of the polylines, i.e., vertices at the ends that are not shared by two edges of a polyline (vertices on the boundary of the structure are ignored). 2. Identify the closest point, from among the neighboring polylines for each isolated point, and record the distance, to this point. 3. Identify isolated points that are within a specified threshold distance of a neighboring polyline, i.e., points with ≤ ℎ ( ℎ was set to 1.5 mm, a value of the same order as fiber spacing, ). 4. Extend the polyline by connecting and if ≤ ℎ . Figure S12a and S12b respectively show the input disconnected fiber layout and the output connected fiber layout generated by this algorithm for the plate with hole structure.

Voxelization Procedure
To build voxel models of the fibrous composites, we developed a specialized voxelization routine in C++ employing OpenVDB along with VTK to support file I/O operations. We chose OpenVDB for its highly efficient sparse dynamic data structure and the associated tools for manipulating these data structures which are designed to efficiently store and work with sparse voxel grids. Voxel grids are essentially 3D images with 3D voxels replacing the 2D pixels. Just as in 2D images, voxels store a value which in OpenVDB's case is determined by the voxel location with respect to the geometry being represented by the voxel grid. The voxels within the volume enclosed by the geometry have negative values and the voxels outside positive.
Our routine takes in as input the geometries of the matrix and any supporting structure in the stl (stereolithography) format and the fiber layout in VTK's vtk or vtp format where the fibers are stored as a collection of polylines. Each of the stl inputs are then directly converted to OpenVDB voxel grids using OpenVDB's in-built function. An empty voxel grid with a background positive value is initially created for the fibers and the voxels are populated with appropriate values iteratively going through each line segment that make up the fibers. The voxels within a distance of (fiber radius) from each line segment are initialized with a negative value, thus creating the voxelized fiber geometry. We perform CSG (Constructive Solid Geometry) Boolean operations as necessary on the voxel grids such that the voxelized matrix, fiber and support geometries have no overlaps. The resulting voxel grids can then be used to directly write slices in bmp format as required by the 3D printer we used. Alternatively, the voxelized model can be written to a file in VTK's meta image format, mhd, for visualization as shown in Fig. S13. In the meta image format, we assign each voxel a unique value to represent the material that occupies the individual voxel. Figure S13a shows the whole 3D voxelized model of a single layer of fibers (magenta) embedded within matrix (light gray), while Fig. S13b shows a slice through the middle of the voxel model.

Material Compiler Validation
To validate the material compilation step, we verified the volume fractions and fiber orientations in the generated fiber layouts for the planar problems. This was achieved with VTK python scripts and the results are shown in Figs. S14-S17. To verify the volume fractions, we randomly sampled 1000 points within the bounding box region of the input geometry. For each point, we calculated the local volume fraction within a 10×10× ℎ mm cuboidal region (ℎ = laminate thickness). The local volume fraction is given by the ratio of volume of the fibers within the cuboidal region to the volume of the cuboid with consideration for points that lie near the boundaries. The calculated volume fractions are shown in Fig. S14a for the plate with a hole problem and Fig. S16a for the laminated plate structure. Figures S14b and S16b show the deviation of the volume fractions in the generated fiber layout from the optimized values in the form of histograms. It can be clearly seen that we get good agreement in both cases. The mean and standard deviation of the error for the laminated problem are 0.07 and 0.76 while they are -0.38 and 1.5 for the plate with a hole structure. The higher deviation for the plate with a hole structure could be ascribed to spatially varying volume fractions as opposed to uniform volume fractions for the laminated plate.
The geometric representation of the generated fiber layout consists of connected line segments, i.e., collection of polylines. The orientation of the line segments thus gives us the local fiber orientation. By calculating the optimal fiber orientation at the segment's mid-point via interpolation, the accuracy of the fiber orientations in the generated fiber layout can be calculated. We defined a quantity, Δ = 1 − |n c ⋅ n o | where n c and n o are unit vectors along the line segment and interpolated optimal fiber orientation respectively. This quantity varies between 0 and 1 with zero signaling no error and 1 implies maximum error (i.e., 90° difference). We calculated Δ for each segment and found that the compiled fiber orientations do not deviate by more than 4° except near defects or singularities as shown in Figs. S15 and S17. It is to be noted that the local material stiffness along the desired fiber direction reduces by 50% when the discrepancy is 4° or more.

Material Characterization Figure S18 Variation of effective modulus, with fiber orientation, of uniaxial tensile testing specimen with the theoretical results obtained from Mori-Tanaka method
To characterize the composite made of the stiff Vero and soft TangoPlus materials, we 3D printed tensile testing samples of length 80 mm and width 20.2 mm on a Stratasys Connex Objet500 multimaterial jetting 3D printer. The samples were 7.7 mm thick such that the volume fraction of continuous Vero fibers embedded in Tango was 10%. We printed seven different samples, each with fibers oriented along a fixed angle with respect to the loading direction which is along the length of the samples. We used the MTS Criterion model 43 universal testing machine with a 1 kN load cell for these tensile tests. The strain applied during the uniaxial tension test was limited to 1% and the strain rate was 1%/min. The effective stiffness modulus along the loading direction for uniaxial tension tests is given by: , where − is the inverse of the homogenized stiffness matrix obtained from the Mori-Tanaka homogenization method. These values for each sample were obtained through a linear fit of the measured stress-strain curves. We used these values to obtain the effective matrix and fiber moduli, and by fitting the experimental data to the Mori-Tanaka method where is a function of and . This was done because we found in our previous studies 3,8 that the directly measured matrix and fiber moduli do not predict the composite properties well, possibly a result of smearing of the matrix-fiber interface during the printing process. Figure S18 shows the experimentally measured values and the best fit ( ) according to Mori-Tanaka method. The best fit obtained by minimizing the difference between experimental values and theoretically predicted values with varying and gave us values of = 0.97 MPa and = 0.85 GPa. These effective moduli were then used in the design automation step of the workflow for the problems where TangoPlus and Vero were considered. Figure S19 Plate with a hole -3D printed optimized designs of different sizes (cases A, B and C from top to bottom) with support tabs to facilitate loading Figures S19 and S20 show the 3D printed planar structures used in our experiments. TangoPlus is known to expand after printing due to strain accumulated during the fabrication process 9,10 . This causes undesired distortion of the printed structure which is visibly apparent in case C given the same thickness but larger size (see Fig. S19).

Digital Manufacturing and Experimental Validation
The optimized structures were printed with relatively rigid tabs (E ~ 1 GPa) to effectively clamp and load the structures. We used the same MTS universal tensile testing machine and 1 kN load cell that we used for material characterization to perform test the printed structures.
For the plate with a hole structure, the tabs were used to clamp in mechanical grips with one of the grips attached to the load cell and the crosshead. The rigid tab also enabled us to impose a uniform displacement load. We specified the rate at which the load was imposed (i.e., crosshead displacement velocity) to be 0.23 mm/min for case A, 0.47 mm/min for case B and 0.7 mm/min for case C such that the average strain rate in the structure remains close to 1 %/min and thus consistent with material characterization. The measured load was found to be linearly proportional and this proportionality constant gives us the effective stiffness of the structure. Given the planar nature of the deformation, we employed digital image correlation (DIC) to obtain the full displacement field from the experiments. The results are shown in Fig. S21 (bottom row) along with the simulation results. The range of displacements from both experiments and simulations agree well.
However, we found the deformation in the y-direction adjacent to the hole to be localized in the simulations but not in the experiments. A high-fidelity simulation of the structure should reveal the source of discrepancy but is beyond the scope of this paper.
The laminated plate structure was simply loaded by placing the printed structure with an elevated supporting frame on a flat plate. A pneumatic gripper attached to the crosshead was used to grip and push the vertical tab downwards at a rate of 0.4 mm/min which equates to an average strain rate of 1%/min in the structure. We used a 100 N load cell given the lower loads needed to bend the plate (~1 N).

Finite Element Analysis with Defects
We set up finite element analysis, while still using the Mori-Tanaka homogenization method, to understand the effect of the singularities/defects generated during the material compilation step. We chose the plate with a hole problem owing to its planar nature and thus low computational complexity. The defects were identified as discontinuities in the fibers as shown in Fig. S15 and located approximately by the open fiber ends (fibers that are too short and fiber ends on the boundaries are ignored). We introduced these defects into the FE model by simply setting the local material volume fraction to zero within a circular region surrounding the open fiber ends while using the optimized fiber volume fractions elsewhere as shown in Fig. S22. The results from these simulations for the three cases A, B, and C are tabulated in Table S3. These FE simulations with defects qualitatively capture the effect of the defects and suggest that the deviation from the ideal performance increases with defect density, thus validating the assumption we made in the main text. When compared to the experiments, we get a good agreement for case B. The reason for the disagreement in case A could be due to insufficient scale separation whereas the disagreement in case C could be attributed to undesired deformations in the printed structure due to excessive swelling of Tango+ matrix material.   Figure S23 illustrates how automated fiber placement (AFP) can, in principle, be adapted for the digital manufacturing step of our workflow. We use the stiff 2D short beam example from Fig. 1 of the main text here. We scaled the structural dimensions by 10 times to 120×80 cm and picked the fiber spacing to be 1 in., a standard dimension for fiber tapes used in AFP and similar processes. Figures S23(a) and (b) show the optimized fiber layouts obtained from the material compilation step with the difference being curves in (b) are offset by 0.5 in. The curves in these layouts can serve as trajectories along which fiber tape can be laid down. Figures S23(c) and (d) show ribbons of width 1 in. obtained by extruding the individual curves along the in-plane normal direction. The ribbons then represent the laid down tapes and it can be seen that we get excellent coverage. The figures also highlight gaps and overlaps occurring near singularities and regions of high curvature. Alternating (a) and (b) layouts during layup will help cover the gaps. High curvature regions could be handled via tow drops. Alternatively, a constraint on maximum curvature could be imposed. This could be achieved by coupling the design automation and material compilation steps and imposing a global constraint on the total curvature.