Design and additive manufacture of functionally graded structures based on digital materials

Voxel-based multimaterial jetting additive manufacturing allows fabrication of digital materials (DMs) at the meso-scale (∼1mm) by controlling the deposition patterns of soft elastomeric and rigid glassy polymers at the voxel-scale (∼90 μm). The digital materials can then be used to create heterogeneous functionally graded material (FGM) structures at the macro-scale (∼10mm) programmed to behave in a predefined manner. This offers huge potential for design and fabrication of novel and complex bespoke mechanical structures. This paper presents a complete design and manufacturing workflow that simultaneously integrates material design, structural design, and product fabrication of FGM structures based on digital materials. This is enabled by a regression analysis of the experimental data on mechanical performance of the DMs i.e., Young’s modulus, tensile strength and elongation at break. This allows us to express the material behavior simply as a function of the microstructural descriptors (in this case, just volume fraction) without having to understand the underlying microstructural mechanics while simultaneously connecting it to the process parameters. Our proposed design and manufacturing approach is then demonstrated and validated in two series of design exercises to devise complex FGM structures. First, we design, computationally predict and experimentally validate the behavior of prescribed designs of FGM tensile structures with different material gradients. Second, we present a design automation approach for optimal FGM structures. The comparison between the simulations and the experiments with the FGM structures shows that the presented design and fabrication workflow based on our modeling approach for DMs at meso-scale can be effectively used to design and predict the performance of FGMs at macro-scale.


U s a g e g u i d e l i n e s
Pl e a s e r ef e r t o t h e u s a g e g ui d eli n e s a t h t t p://s u r e . s u n d e rl a n d. a c. u k/ p olici e s. h t ml o r al t e r n a tiv ely c o n t a c t s u r e @ s u n d e rl a n d. a c. u k.

Introduction
Additive manufacturing (AM) based on homogeneous metal, ceramic or polymer materials has made inroads into manufacturing environments [1] such as aerospace, medical devices, spare-part applications, and short-run production environments [2,3]. These industries take advantage of the increased geometrical complexity [4], and the mass-customization capabilities [5] of AM. However, the true advantage of AM technologies in product design and manufacturing is when: "Product performance is maximized through the synthesis of shapes, sizes, hierarchical structures and material compositions, subject to the capabilities of AM technologies" [6].
In this regard, selective deposition of heterogeneous materials offers unprecedented possibilities by enabling the fabrication of 3-dimensional (3D) objects composed of multiple materials with different physical and/or chemical properties [7]. Digital fabrication of complex structures composed of functionally graded materials (FGMs) is now possible [8]. FGMs display complex spatially varying material properties [9]. For example, combinations of rigid and soft materials can form custom-designed anisotropic properties or property gradients that cannot be generated otherwise within a single build using a single material system [10]. Controlled deposition of multiple materials provides extensive capabilities in engineering design by enabling the customization of material microstructures and therefore, programming of the part to achieve prescribed functionality [11][12][13][14][15][16][17].
Multimaterial AM systems with the ability to vary material compositions locally during the build show great potential for FGM applications. AM technologies are evolving in this direction, examples include vat-photopolymerization, material extrusion, binder jetting, directed energy deposition, powder bed fusion, and material jetting process categories [18,19]. The design and manufacturability of AMderived FGMs is limited to how discretely the AM system can selectively deposit the base materials [20]. Among the commercially available multimaterial AM systems, material jetting 3D printers are known for their high resolution [21,22].
Material jetting incorporates the widest variety of colors and materials into a single manufacturing system among all AM technologies [20]. Due to the discrete and digital nature of the deposition system, the multimaterial system obtained through this process is termed digital materials (DMs) [23]. DMs, in our context, are effectively homogenous materials made of discrete elements of two or more isotropic materials. The basic principle of the material jetting technology is to use a block of piezoelectric inkjet heads with multiple nozzles that are used to deposit droplets of multiple photopolymer resins following computer-generated deposition patterns. The droplets are cured and flattened to a sheet by an ultraviolet light source and a roller respectively, both adjacent to the inkjet heads. These cured and flattened droplets form the basic building blocks i.e. volume elements or voxels which when put together can realize any solid freeform subject to the usual restrictions of support material removal. Material jetting 3D printers, such as Stratasys J750 which we used in this study, allow production of multimaterial 3D structures using user-defined voxel patterns (called voxel printing) in addition to the standard practice of printing with STL (stereolithography) files.
The mechanical behavior of parts made by material jetting can be defined by the combined effect of mechanical properties of the base materials used at the voxel-scale (∼90 μm). This multimaterial system once realized can then be viewed as a new composite or a so-called digital material at the meso-scale (∼1 mm) with properties different from the constituent materials at the voxel-scale. Mathematically, there is no unique arrangement or pattern of voxels to obtain a desired material response. Stratasys provides pre-defined proprietary patterns for a few digital materials through the standard printing mode [21]. However, one can create a continuous gamut of digital materials with the help of voxel printing [10,24,25]. This gives extensive design freedom to create DM-based FGMs as one can now design digital materials with desired properties and gradations within the physical limits imposed by the choice of the constituent base materials. The design of FGM structures requires the simultaneous integration of material design, structural design and product fabrication [11,22].
Among the design automation based FGM design approaches, some [30,31,34,36] do not connect the FGMs to any specific microstructure while others [32,33,37,38] connect them to lattice or cellular microstructures. Weeger et al. [24], and Zhu et al. [25], are two of the design automation efforts that utilized DM-based FGMs. In any FGM design approach, the prediction of the overall mechanical response at macro-scale requires connecting the meso-scale design variables such as volume fractions and orientations to the macro-scale properties such as Young's modulus and density either through concurrent numerical modeling of the micro-scale or high-fidelity material models, typically based on mathematical homogenization.
Mathematically, there is no unique way to connect material gradation to the underlying DMs. Hiller and Lipson [27], in their work demonstrate different ways of doing so including random, mesh, layer, dither and longitudinal material distribution strategies. To obtain isotropic material properties, random distribution [24,45] and dithering [10,29,46,47] approaches were used. Zhu et al. [25] used inverse homogenization while Weeger et al. [24] used random dithering for DM-based FGM design.
Most studies on the characterization of DMs use the standard printing mode of material jetting machines and do not explore the full capabilities of voxel printing. The standard printing mode refers to material combinations provided by Stratasys where the volume fractions of soft and rigid polymers and the deposition patterns are hidden from users. These material combinations are preconfigured in the preprocessing software [20]. Previous research has used this standard printing mode to study the mechanical performance of DMs in terms of anisotropic behavior of the elastic modulus, the ultimate tensile strength, and the percentage of elongation at break with varying process and material ageing conditions, specifically for rigid Vero-WhitePlus resins in [48,49] and combinations of DM in [50]. Salcedo et al. [8] fabricated and tested tensile specimens with embedded regions of FGMs comprised of the pre-defined DMs to generate and validate material models. Similarly, different hyperelastic material models were used to simulate the tensile behavior of predefined DM compositions by Ryu et al [51]. It was found that the material characterization of the composite DMs exhibit significant non-linearity under large strains with a rate-dependent behavior [52], with properties similar to conventional short-fiber composites [53]. Altogether, current material models are limited in their ability to describe the manufacturing process and microstructural dependency to accurately predict the behavior of DMs and thus affect the design of FGM components. A need for a probabilistic approach to material modeling that accounts for uncertainties was proposed [54]. Such an approach could aid in effective design of FGM structures by considering the processstructure-property relationships and the underlying uncertainties.
Designing and simulating FGM parts is challenging, especially when the properties of the graded material vary rapidly in space. In this case, grading is reflected at both meso-scale (∼1 mm) and macro-scale (∼10 mm) i.e., the overall structural scale [54]. To design FGMs, we need to effectively describe the microstructural behavior of DMs with a representative material model, which can be used to effectively simulate the material behavior arising from the non-uniform material distributions. The standard approach based on homogenization procedures entails identifying various microstructural descriptors (such as volume fraction, geometry and topology of the arrangement of the constituents in the DMs) and understanding the underlying mechanisms that tie the microscale mechanics to the macro-scale structural behavior. Homogenization procedures generally preclude process-induced uncertainties or variations such as anisotropy or defects.
We adapted here an approach based on regression analysis where we fit the DM characterization results to a cubic regression model to predict the macro-scale behavior of DM-based FGMs. The use of regression analysis carries an inherent advantage of not only capturing the microstructural behavior but also the process-structure-property relationships. In our model, the microstructural descriptor is the volume fraction of the rigid material within the DM and the process parameter is printing orientation which accounts for process-induced anisotropy. Regression models could also lead us towards a probabilistic approach to material modeling owing to the ability to capture prediction intervals representing a range of likely values or variation for new observations.
In summary, this paper presents a design and fabrication workflow for DM-based FGM structures. To facilitate the material modeling and thereby the design process, we fabricated and tested DM tensile specimens each with different volume fractions (microstructure descriptor) and print orientations (process parameter) through random yet controlled deposition patterns of rigid polymer, VeroClear (E˜1 GPa), and soft elastomer, TangoPlus (E˜1 MPa) [15]. Based on the material characterization results, we created a series of regression equations to describe the deformation and failure properties of the DMs including the effect of the printing orientation to account for process-induced anisotropy. We used the regression-based material model to then demonstrate with simple planar structures a workflow that integrates design and fabrication of DM-based FGM tensile structures with varying material gradients. We employed two different design approaches. In each case, the designed structures were fabricated and tested to validate the material model and the workflow. In effect, we present here a general methodology for characterization, modeling, design and fabrication of DM-based FGM structures making it one of the first few efforts to connect material characterization and FGM structural design to fabrication and validation.

FGM design and manufacturing workflow
In Fig. 1, we present the developed workflow with its four key steps identified. The workflow begins with the engineering problem specification i.e., structural geometry, load, boundary conditions and design goals along with the relevant material model for DMs. This material model could be based on theory, simulations or experiments. Here, we relied on an experiment-driven regression-based model that connects mechanical properties to microstructural descriptors and potentially, process parameters. The DM characterization experimental setup and the results are presented in Section 3.
The second step of the workflow deals with FGM structural design wherein the design goal is achieved by either: a) prescribed (see Section 4.1) or b) automated design approach (see Section 4.2). In the prescribed design approach, we used intuitive gradation patterns for the material properties and applied them to structures to tailor their mechanical response. For the automated design approach, we employed a TO method to generate optimal FGM structures. The FGM realization step helps translate DM-based FGM designs to manufacturable distributions of rigid and soft polymers, thus connecting the design step to the final step of product fabrication via voxel-based multimaterial jetting (see Section 3.1 for details). The workflow enables seamless integration of design and fabrication while also enabling simultaneous design of material at meso-scale through DMs and structure via gradation of DMs.

Digital materials fabrication and process planning
The material jetting AM system used for the experiments is Stratasys J750. The machine consists of a block of eight inkjet print-heads that feature full-color printing capabilities and combine up to five materials in the same printing process. The material jetting process consists of inkjet-based deposition of low viscosity photopolymer resins onto the build plate, creating an array of droplets and the object layer by layer. Each layer is cured by an ultraviolet lamp directly after deposition and then flattened with a mechanical roller in a continuous process. Each inkjet head is composed of 96 circular nozzles arranged in a linear array, each measuring 0.5 mm in diameter [55]. The two photopolymers used in this work were the flexible TangoPlus (FLX980) and the rigid VeroClear (RGD810). For simplicity, VeroClear and TangoPlus will be referred to as Vero (V) and Tango (T) respectively. We defined the fraction of Vero by volume in a unit volume of the DM as volume fraction, f .
To produce DMs with different volume fractions, we used the voxel printing capability which gives us the ability to program the material jetting process using a set of pixelated images in PNG (Portable Network Graphics) format. Each image refers to a horizontal slice of the 3D model to be printed that is parallel to the build platform and has one layer of voxels in the normal direction. The image's pixels refer to locations of voxels within this layer and are each assigned an arbitrary color. These colors are mapped before 3D printing to the specific materials to be jetted with one of the colors representing void. The image resolution and the number of slices are commensurate with the bounding box of the 3D printed model and the voxel size. Thus, a DM with a desired volume fraction, f is obtained by randomly assigning the voxels in each slice to Vero (Tango) with a probability of f (1f ). Fig. 2 shows the schematic representation of the key machine elements as well as the principles of the voxel-based material jetting process: (a) represents the printing block which moves in the X and Y directions to deposit the resin, (b) is the build platform which is responsible for the motion in Z direction, (c) is the ultraviolet lamp responsible for curing the photopolymer resin, (d) shows the 8 printheads of the print-block, (e) shows the representation of the roller that is used to flatten the deposited layer of photopolymer resin, (f) represents the DM specimen under fabrication, (g) shows a section of the pixelated image corresponding to a slice of the 3D model and a zoom-in depicting the voxels and (h) is the layout of samples with different build orientations.
In Fig. 2(g), the white and gray pixels correspond to Vero and Tango respectively. In addition, a schematic representation of the voxels is also shown. The smallest voxel size achievable with Stratasys J750 is 42.33, 84.66, and 13.5 μm in the X, Y and Z directions, respectively [56]. A custom voxel size could be chosen but the voxel dimensions in X, Y or Z directions must be an integer multiple of the minimum voxel sizes. We chose a recommended voxel size of 42.33 (V x ) x 84.66 (V y ) x 27 (V z ) μm for faster print times while retaining the high print resolution.
As part of the FGM realization step of the workflow (step 3 in Fig. 1), an extended cellular automaton-based pseudorandom number generator as implemented in Mathematica was used to create the image slices that comprise the DM-based FGM structures. As done with DMs, each pixel was associated with Vero (Tango) with a probability equal to the local volume fraction, f ( 1-f ). This method allows us to distribute randomly yet in a homogeneous fashion the two photosensitive resins while guaranteeing a prescribed volume fraction at the meso-scale (∼1 mm). We can thus design FGMs at the macro-scale (∼10 mm) with varying material properties by distributing the DMs according to a prescribed distribution or gradation of volume fractions, f x y z ( , , ).

Digital materials characterization
We prepared and characterized, via voxel printing, several rectangular DM tensile specimens (see Fig. 2(h)), each with a uniform volume fraction to obtain their mechanical properties. Specifically, we employed a design of experiments (DOE) method to measure Young's modulus, elongation at break and ultimate tensile strength and determine their dependence on the microstructural descriptor (i.e. the volume fraction) and the manufacturing process parameter (i.e. the printing orientation).
We tested six Vero volume fractions ranging from = f 1, which represents 100% Vero to = f 0, which represents 0% Vero or equivalently 100% Tango. In addition, the samples were oriented in three different build orientations XY, YX and ZY as shown in Fig. 2. Both XY and YX are parallel to the build plate, whereas ZY is oriented in the perpendicular direction. The dimensions for the tensile specimen printed horizontally (i.e. XY and YX) were 25 × 5 x 60 mm with a cross-sectional area of 125 mm 2 , whereas the specimens oriented vertically (i.e. ZY) were 25 × 5 x 50 mm with the same cross-sectional area of 125 mm 2 . In the case of the ZY tensile specimens, the gauge length was reduced by 10 mm to limit the printing time. The DM characterization was conducted using a uniaxial tensile test on an Instron model 5982 mechanical testing machine with a 100 kN load cell. With a pneumatic gripping force of approximately 10 kN, the specimens were elongated until failure at a strain rate of 0.4 mm/sec. These characterization data were used to build the regression models for Young's modulus (E), ultimate tensile strength (UTS), and elongation at break (EAB) as a function of f per each build orientation XY, XZ, and YZ. Table 1 shows the two independent variables used to create the regression model for the mechanical properties of the multimaterial Vero-Tango system.
As part of the DOE, an Analysis of Variance (ANOVA) test was conducted. The ANOVA of the collected data was carried out at a confidence level of 95% (α = 0.05). The ANOVA test included both independent variables: f and sample orientation (Orient.). During the analysis, f was treated as a continuous factor subject to 0 f 1; whereas, sample orientation was treated as a categorical factor, ultimately creating one equation per sample orientation XY, YX, and ZY. Eq. (1) defines the scalar form of the implemented third order polynomial regression model: Here, is the unobserved random error, with appropriate subscript represents the coefficients of the regression model for each term (i.e. intercept , first order i , quadratic ii , and cubic iii ). These are calculated by the least squares method. x i represents the independent variable volume fraction f , x j represents sample orientation, which is treated as a categorical variable, and y corresponds to the dependent variables (i.e. Young's modulus, ultimate tensile strength, and elongation at break).
To construct the regression models, all first, second, and third order terms were included in the test. The statistical significance of the terms was assessed by looking at the test's p-value. If the p-value was below the specified significance level (i.e. α = 0.05), the term was declared to be statistically significant and the test's null hypothesis was rejected [57]. These best fit models were later used to model the DMs in the FGM structural design step of the workflow (step 2 in Fig. 1).

Digital material model
We begin with a brief discussion of the uniaxial tensile test results before presenting the regression-based material model for DMs. Fig. 3 shows the stress-strain curves obtained during the uniaxial tension tests of the DM tensile specimens with different volume fractions and oriented along XY. The results show that the specimens with higher volume fractions are highly nonlinear while those with low f behave linearly with an increased ability to elongate up to 130% of the initial length.
One may see that the ultimate strain decreases with the increase of Vero volume fraction. On the contrary, ultimate tensile strength   Fig. 3. The stress-strain curves of the XY oriented tensile specimens with a volume fraction: f 0 1 subject to a strain rate of 0.4 mm/sec. Also shown are the representative slices of the respective tensile specimens to illustrate the composition of the tested DMs.
increased with increased Vero volume fraction. For specimens with f ≥ 0.8, the results show higher ultimate tensile strength at the cost of a brittle fracture and limited capability to elongate. Specimens within the range of 0.4 f 0.6 show a strain hardening phenomenon, where the DM material is strengthening by plastic deformation after the elastic limit is reached. Finally, the specimens with f 0.2 show a linear elastic deformation until the moment of fracture. This phenomenon is consistently reproduced independently of the build orientation. However, the measured ultimate tensile strength (UTS), elongation at break (EAB) as well as Young's modulus (E) are orientation dependent and consistent with the existing research [50]. We refer to Appendix A for a detailed discussion regarding our findings on the anisotropic behavior of the DMs.
As mentioned in Section 3.2, we performed the ANOVA test treating f as a continuous factor, with a range of values: f 0 1 and printing orientation as a categorical factor for XY, YX and ZY. The results of the ANOVA (more details in Appendix B) showed that all three first, second, and third order terms were statistically significant; and therefore, they were included as terms in the regression models. Table 2 shows the model summary in terms of standard deviation (S) that accounts for the distance between the data values and the fitted values and percentage of variation (R 2 ) which is used to determine how well the model fits the data. The cubic regression models in natural units for each dependent variable as a function of f for all three orientations XY, YX, and ZY are also presented in Table 2.
An example fitted curve corresponding to the Young's modulus (E) as a function of the Vero volume fraction for specimens oriented along XY is displayed in Fig. 4. The figure also displays the confidence interval (CI) at 95% that represents a range of likely values for the mean response as well as the prediction interval (PI) at 95% that is the range of likely values for new observations.
It should be mentioned that the presented regression model fails to calculate the real Young's modulus of the DMs with volume fractions f ≤ 0.2. In this region, the regression model would provide a negative Young's modulus. Ideally, it is possible to use a piecewise function with two regions that define E f orient ( , .) for 0 f 0.2 and 0.2 <f 1.0 respectively, thus transforming the data into a regression with reduced lack of fit, and assuring that the prediction of E f orient ( , .) has positive values. To avoid numerical issues and to keep the regression models simple, we limited ourselves to f ≥ 0.2 when using the model in the design of FGM components.

Design of FGM structures
Here, we demonstrated and validated our FGM structure design and manufacture workflow discussed in section 2. For the design problems, we limited ourselves to linear elastic planar structures and employed the regression model of Young's modulus along the XY orientation presented in Fig. 4. We used two design approaches to engineer FGMs, step 3 of our workflow. In the first approach termed prescribed design approach, we designed uniaxially-loaded tensile structures by simply prescribing various graded DM distributions.
In the second approach that we called automated design approach, we optimally designed FGM structures via the design automation method of topology optimization. The regression-based material model was embedded into a topology optimization problem with the goal of obtaining an FGM structure that achieves a predefined displacement map and thereby a desired mechanical response. To demonstrate the versatility of the approach, two different problems were formulated with the same objective but with different loading conditions. The FGM designs from both approaches were realized and printed, thus exercising the full workflow. The workflow and the material model were validated by comparing the FE analysis predictions with those from experiments.  Fig. 6 shows the representative slices of voxel patterns generated to realize the material gradations shown in Fig. 5 as part of the FGM  realization step of our workflow. This is achieved by placing Vero at a voxel location, (x , y) with a probability dictated by the chosen gradation function f x y 0 ( , ) 1. This creates locally a digital mixture of Vero and Tango materials consistent with the given gradation patterns as shown in Fig. 5. This is done for each slice along the thickness direction of the tensile FGM structure. The effective modulus of each structure was experimentally measured following the same procedure detailed in Section 3.2. The effective modulus, here, refers to the force per unit area required to support the specimen under unit tensile strain.

Prescribed design approach
The FE simulation results were compared with the experiments to validate our approach. We assumed a constant material Poisson's ratio of 0.388 [54] independent of f and print orientation. We used plane stress 4-noded quadrilateral elements of size 0.5 mm to model the structure and imposed clamped boundary conditions on one of the short edges while imposing a unit displacement boundary condition on the other to replicate the experimental conditions. The material properties in each element were varied in accordance with the gradation function, f x y ( , ). Specifically, the Young's modulus at a location (x , y) was varied according to: The effective modulus was then calculated by obtaining the total reaction force at the clamped edge per unit area and dividing it by the strain induced by the unit displacement (i.e. = x 0.01667). Fig. 7 and Table 3 show the experimental results and comparison with FE simulations of graded tensile structures respectively. Fig. 7 shows the stressstrain curves for each case (3 curves for each of the 3 printed specimens). All the structures exhibit a stiffer response at small strains (< 0.1%) and a softer response thereafter. We fit the response to a linear function, the slope of which gives us the effective modulus of the specimen as tabulated in Table 3. The measured values compare well (within 15.05%) with the FE results in all the four different cases as can be seen in both Fig. 7 and Table 3, thus validating the regression models and our workflow with the prescribed design approach. It is to be noted that in the case (c) and (d) in Fig. 6, the gradation patterns lead to coupling between tension and bending deformations. In the experiments, we only measured the longitudinal forces and ignored the bending forces while the coupling is implicitly accounted for in the FE simulations.

Automated design approach
To further validate the material modeling approach and demonstrate our workflow with the ability to automate the design of optimal FGM structures, we conceived a design automation problem based on optimization with the goal of obtaining an FGM structure that achieves a prescribed longitudinal displacement map i.e., = u x y u x y ( , ) ( , ) Here, s are the optimization variables used to define the local Vero volume fraction, = f x y f x y ( , ) ( , ) (i.e., designs are symmetric about xaxis). This was done in terms of a standard linear filter with a filter radius of 2 mm to regularize the optimization problem. V vero and V are the total volume of Vero within the design domain, , and the total domain volume respectively while f is the Vero volume fraction in the design domain. The objective, z, is a function of the x -displacements, u x , which in turn are functions of the spatial coordinates, x and y, and the optimization variables, s. This relationship is a result of the underlying physics i.e. elastostatics described by the equation: = u ( ) 0, where is the 4th-order stress tensor defined = C E f s u ( ( ( )), ) ( ) with C being the isotropic material stiffness tensor and , linear elastic strain tensor.
The elastic modulus was defined in terms of f in accordance with our experimental findings (see Eq. (2)). We solved two sets of optimization problems with an identical setup except for the loading conditions. In each case, the objective was to minimize the squared error between the prescribed displacement map in the x-direction, = u u x l ( / ) . This guaranteed that as the objective value was reduced through the course of the optimization procedure, we obtained a design whose longitudinal deformation behavior got progressively  x p , with = z 0 being an exact match.
In the first set of optimization problems, hereafter referred to as "problem 1", the load was a prescribed uniform displacement, = u 1 0 mm in the x-direction, while the y-displacements were set to zero along the right edge of the structure, see Fig. 8(a). In the other set of problems, hereafter referred to as "problem 2", the prescribed displacements were limited to the highlighted portion of the right edge as shown in Fig. 8(b). The other end was rigidly clamped i.e. = = u u 0 x y in both problems. With an isotropic material, such a loading condition results in a nearly uniform longitudinal strain distribution or in other words a linear variation of longitudinal displacements, u x . Our target, however, was to obtain a linear variation of the longitudinal strain i.e., a quadratically varying displacement in the x-direction as shown in Fig. 8(c). This can only be realized by a non-uniform distribution of DMs and thereby Vero and Tango. We solved the problems using topology optimization where the objective is evaluated through FE analysis, a search direction is arrived at through sensitivity analysis and the design updated by an optimization algorithm, GCMMA (globally convergent method of moving asymptotes) [58] in our case. This constraint on f effectively lets us dictate the stiffness of the structure as higher Vero content leads to stiffer structures while also achieving the principal objective of a linearly varying longitudinal strain. Table 4 outlines the optimization results where the final optimal objective values, z and the effective stiffness of the structure (defined as the ratio of reaction force to displacement, u o ) are presented for each case of problems 1 and 2. We get consistently better results with problem 1 for each case when compared to the equivalent case of problem 2. This is to be expected because in problem 1, the error on the right edge is zero by design whereas in problem 2, this is not so for the free portion of the right edge. This results in significant material gradation in both x-and y-directions in problem 2 but not observed in problem 1. For problem 2, we observed the creation of a stiff mechanism that keeps the x-displacements nearly rigid in order to achieve the target displacements locally. We get the best results in each problem with the unconstrained case where the final Vero usage is 28% for problem 1 while it is 38% for problem 2. The case with = f 0.8 performs the worst due to the excessive use of Vero. In general, we get a higher stiffness in each of the constrained cases at the expense of the principal objective.
We manufactured the designs obtained from the unconstrained cases using the same method as the one used for prescribed FGM designs where the local volume fraction at a material point was used to define the deposition patterns. For the mechanical testing of these FGMs, we used an MTS Criterion Model 43, equipped with mechanical grippers, 10 kN load cell and a digital image correlation (DIC) system. The DIC system was used to capture the displacement maps of the tested specimens. We used an Allied Vision, Manta G-146 vision camera with a maximum frame rate of 17.8 fps at full resolution and a general-purpose lens with a focal length of 8 mm. The specimens were sprayed to obtain an optimal speckle pattern and evenly illuminated to avoid overexposed highlights on the target surface, thus allowing precise measurement of displacements. Fig. 11 shows the comparison of results obtained from the FE  simulations and experiments with DIC. Fig. 11(a) and (b) show the representative slices of the realized optimal FGM design (magenta -Vero, white -Tango) obtained from solving the optimization problems 1 and 2 respectively while Fig. 11(c) and (d) show 3D printed translucent mockups of the optimal structures, printed with transparent and magenta colored materials respectively in place of Tango and Vero and with volume fractions scaled down by 10 times to achieve the translucency. Fig. 11(e) and (f) show the longitudinal displacement maps obtained from FE simulations for problems 1 and 2, respectively. For problem 1, we get good agreement with the prescribed displacement map (see Fig. 8(c)), whereas for problem 2, the obtained displacement map deviates considerably near the loaded right edge while agreeing well elsewhere. Fig. 11(g) and (h) show the experimentally obtained displacement maps via DIC. We get excellent agreement pointwise between the DIC experiments and FE simulations thus validating our DMbased FGM design approach. We also measured the effective stiffness of the printed structures along with the displacement maps. We obtained values of 241.22 and 416.39 N/mm for problems 1 and 2 respectively. These values compare well with the simulations where we get values of 279.28 and 410.12 N/mm for the same and offer further validation of our approach.

Conclusions
We presented a design to fabrication workflow for DM-based FGM structures that integrates material and structural design while digitally connecting them to fabrication. We demonstrated and experimentally validated the workflow with simple planar structures using two different approaches: prescribed and automated. In the prescribed design approach, designer specified gradation patterns were directly applied to FGM structures while in the automated approach, optimal gradation patterns were obtained via a design automation technique. To facilitate efficient design of the FGMs, we made use of regression analysis of the experimental results obtained from the mechanical characterization of DMs. The use of regression analysis not only simplified the prediction of the effective mechanics of the DMs using the microstructural descriptors ( f ) but also allowed us to account for the process related parameters (orient.) that affect the structural behavior.
In the case of the prescribed FGM designs, the maximum absolute error between the FE simulations and the experiments was 15.05% while it was 15.78% in the case of automated FGM designs. The sources of these discrepancies include uncertainties arising or introduced during fabrication, material characterization, and experimental validation steps. Future research employing probabilistic material models that can propagate the uncertainty to FE simulations can mitigate the discrepancies observed here. Nevertheless, the comparison between the displacement maps obtained via FE simulations and DIC experiments of the automated designs showed an excellent pointwise agreement.
In summary, this research is a first step towards development of integrated design to fabrication workflows for FGM structures based on multimaterial AM systems. Future work that accounts for anisotropy, non-linearity and failure will allow us to describe the elastic behavior of DMs with better accuracy and enable us to engineer FGMs useful for 3D structures such as energy absorbers, soft robotic actuators or similar complex devices that can benefit from locally varying the stiffness behavior to achieve the desired performance.
While the focus of this paper was restricted to mechanics and an objective of matching the displacement field to a target, this approach can be extended to more diverse and sophisticated objectives as dictated by a particular application e.g. damping and energy absorption mechanisms, tough fracture resistant structures, thermo-mechanics, acoustics, shape memory materials, stimuli-responsive materials, and 4D printing.
Given this potential, future research will ideally extend this method   to other multimaterial AM process categories, such as (i) directed energy deposition with multiple metal alloys, (ii) powder bed fusion for metals or ceramics through porosity, micro-structure customization and combination of materials in the build process, (iii) multimaterial nanoparticles material jetting, (iv) material extrusion of FGM through path planning and dual filament extrusion, as well as (v) binder jetting through binder concentrations.

Declaration of Competing Interest
None.

Funding
None.
To account for the anisotropic behavior of the DMs, specimens were positioned on the build platform in three different orientations as introduced in Fig. 2 and Table 1. Fig. A1(a) shows the stress-strain curves for XY, YX and ZY oriented specimens with = f 1 and Fig. A1(b) shows the same stressstrain curves for the specimens with = f 0. In general, for increased f the material behaves like a rigid polymer with a higher degree of anisotropic behavior. The results for higher f show that XY and YX oriented specimens have higher UTS and EAB in comparison to specimen oriented along ZY. At the same time, E varies with the build orientation. Parts oriented vertically in the Z-axis are the weakest. On the opposite end, for decreased f the material behaves elastically and is isotropic in terms of stiffness E. However, UTS and EAB results are anisotropic and therefore orientation dependent.
Figs. A2-A4 show the effect of increasing f on Young's modulus (E), ultimate tensile strength (UTS) and elongation at break (EAB), respectively. The three figures also display the results for XY, YX and ZY orientations, and each data point represents the mean value of three repetitions per Fig. 11. Optimization, Fabrication, FE simulation and DIC experimental results. Representative slices of the realized optimal FGM design (magenta -Vero, white -Tango) obtained from solving the optimization problems 1 and 2 are respectively shown in (a) and (b) while in (c) and (d), 3D printed translucent mockups with clear and magenta colored materials are shown for visualization. Maps of x -displacements obtained from FE simulations for problems 1 and 2 are shown in (e) and (f) respectively whereas (g) and (h) show the results obtained from the DIC experiments respectively. experimental combination and the standard error. Evidently, the standard error for the measurements is found to be negligible. Fig. A2 shows that increasing Vero volume fraction, f results in increasing Young's Modulus (E). Likewise, the anisotropic effect induced by the printing orientation increases with f . For specimens oriented along XY, YX and ZY, the range of E (MPa) is between 0.53 E XY (Orient.) 1401.67, 0.56 E YX (Orient.) 1206, and 0.46 E ZY (Orient.) 1083.33 respectively. The relative difference between the weakest and strongest printing orientations is 13.2% and 22.7% for the soft Tango ( = f 0) and rigid Vero ( = f 1) respectively. The results in Fig. A3 show that orientation has significant impact on UTS. With increasing Vero volume fraction, UTS as well as its anisotropy increase. For specimens oriented along XY, YX and ZY, UTS (MPa) ranges between 0.51 UTS XY 60.33, 0.53 UTS YX 59.72, and 0.40 UTS ZY 38.66 respectively. The relative difference between the weakest and strongest orientations is 21.6% and 35.9% for the soft elastomer ( = f 0) and the rigid polymer ( = f 1) respectively. The effect of orientation is more significant for EAB, Fig. A4 show how specimens oriented in the weakest direction ZY have a decreased ability to elongate in comparison to XY and YX orientations. In average terms, the relative difference for EAB between the weakest and strongest build orientations is 57.6%. In general, the anisotropic behavior is equally significant for specimens with low or high f and the EAB decreases as f increases. In summary, the results show that all the measured mechanical properties are found to be weak with specimens manufactured vertically in the ZY orientation due to the layered nature of the AM process. The strongest mechanical properties are found to be in the XY orientation, which is perpendicular to the print-heads of the print-block as shown in Fig. 2. Each DOE combination presented in Table 1 was repeated 3 times with a total of 54 tensile test runs. Table A1 shows the averaged measurements for Young's modulus (E), ultimate tensile strength (UTS) and elongation at break (EAB) with different Vero volume fractions and print orientations. Table B1 shows the full ANOVA table for first order, second order, and third order. Additionally, bilinear and quadratic interaction terms are also included in the test. To understand interaction and significance of the terms, we can compare the variance of the effects to the variance of the residuals. We applied a backwards elimination algorithm of non-significant terms that removes the least significant terms at each step. The algorithm stops when all the terms in the model have p-values that are lower or equal to 0.05. This provides the best fit model to represent mechanical properties of DMs in terms of the microstructural descriptor, volume fraction and the manufacturing process parameter, the printing orientation.

Appendix B
The ANOVA table provides this information in the form degrees of freedom (DF) which accounts for the amount of information of the data set in terms of replicates for each observation. The results show a much higher DF for the pure error term of the ANOVA table, and therefore, a good fit of the regression model to predict the response of the dependent variables. Additionally, the F-value is sufficiently large indicating that the first, second, and third order terms of the model are significant, and need be included in the regression model.
In relation to the statistical significance of the orientation and volume fraction f ( ), the p-value is used to show the probability that measures the evidence against the null hypothesis. The lower the p-value, the stronger the evidence against the null hypothesis. In this case, the null hypothesis is explained as the lack of significant difference between specified populations; any observed difference is due to sampling or experimental error. As a result, first, second, and third order terms related to f show great statistical significance in the variation of the dependent variables E, UTS, and EAB with p-value ≤ 0.1; thereby, rejecting the null hypothesis. Furthermore, the ANOVA test shows that the effect of orientation is statistically significant for the variation of all three dependent variables. Consequently, showing p-values of 0.766, 0.433, and 0 for E, UTS and EAB, respectively.
The effect of orientation is especially significant for EAB when compared to UTS and E. The ANOVA table shows that bilinear and quadratic interaction terms are significant for EAB and UTS. On the contrary, only the bilinear interaction term between f and orientation is significant for E. Ultimately, the regression equations presented in Table 2 do not integrate interaction terms as the orientation was defined as a categorical variable with three values XY, YX, and ZY. This is a first step towards characterization of DMs to 3D design problems that can benefit from FGM use. Future characterization is necessary to study the anisotropic elasticity of shear modulus and Poisson's ratio as a function of Vero volume fraction f and orientation. This work will allow description of the stiffness matrix of DMs using an orthotropic linear elastic material model.