Computational Design of Compositionally Graded Alloys for Property Monotonicity

Materials that leverage spatial property gradients or functionally graded materials (FGMs) have been investigated for decades [1,2], but are now a rapidly expanding area of research due to the emergence of additive manufacturing [3–5]. These materials are typically used to achieve contradictory performance requirements in the same part. For example, metal-ceramic composites trade mechanical strength for thermal resistance [1], while graded polymers balance density and flexibility [5]. FGMs can be created by varying material composition, processing, or structure and have been demonstrated in polymers, ceramics, and metals [3–5].

Compositionally graded alloys are a special class of FGMs that exhibit the high performance characteristics of metals, but offer larger and more varied changes in properties than grading structure alone [6–9]. Directed energy deposition (DED) [10] greatly simplified the construction of such materials by enabling manufacturers to change the deposited composition layer by layer. During the DED process, powders of various compositions are deposited on the build and then solidified by a high energy laser. The composition of the build can be easily controlled by controlling the ratio of deposited powders.

In recent years, many compositionally graded alloys have been printed using the DED process [11–16]. While some of these works have been successful, almost all linearly grade composition directly between two material endpoints. As such, they often encounter deleterious phases in the gradient path that can lead to undesirable properties or cracking during the build process. Carroll et al. [12] and Chen et al. [13] attempted linear gradients between 304L stainless steel and Inconel 625, but observed secondary phases that led to cracks and increased microhardness. Reichardt et al. [14] attempted a gradient from 304L stainless steel to Ti-6Al-4V with a pure V interlayer, but experienced significant cracking due to the formation of the brittle sigma phase. Bobbio et al. [15] linearly graded between Ti-6Al-4V and Invar and also saw cracking due to the formation of brittle intermetallics. Meng et al. [16] manufactured a linear gradient from Ti-6Al-4V to Inconel 625 that, once again, cracked due to the formation of secondary phases.

Hofmann et al. [6] proposed a strategy to avoid the formation of these detrimental secondary phases. In that work, the authors proposed using phase diagrams as maps to plan gradient paths so that undesirable phases are avoided. The strategy was employed in Reichardt et al. [14] when CALculation of PHAse Diagrams (calphad) software was used to generate an isothermal ternary phase diagram in the Fe-Cr-V system. This phase diagram was then used to plan a new gradient from pure V to 304L stainless steel that avoids the sigma phase region the authors encountered in the linear gradient. Recent works have built on this strategy by using Scheil ternary projections, which provide more accurate predictions of undesirable phase regions that consider the rapid solidification conditions of additive manufacturing [17,18].

While this strategy is useful for visualizing and planning simple gradients or sections of more complicated gradients, it has several significant limitations. The biggest impairment is the inability to visualize more than three elements at a time when most gradient design problems involve many more elements. For example, a gradient between 304L stainless steel and Ti-6Al-4V will necessarily contain significant quantities of at least six elements (Fe, Ni, Cr, Ti, Al, and V) which have 20 possible ternary combinations. It would be impossible to sufficiently visualize this entire composition space with any two-dimensional illustration. Such phase diagrams are also isothermal or temperature projections, limiting the ability to account for all the temperatures a gradient might experience during manufacturing or operation. Lastly, phase diagrams do not provide information about the property space, so designers have no ability to ensure desirable properties of a designed gradient.

In Kirk et al. [19], a new computational FGM design methodology was presented that encodes some of these design principles computationally. In short, a surrogate model of phase regions is trained from calphad software and then a path planning algorithm is used to design paths that avoid undesirable phase regions while optimizing a specified cost function. As a computational methodology, it escapes the limitations of human visualization and enables the design of gradient paths in high dimensional spaces. It has been shown to avoid undesirable phase regions and even optimize gradient paths for cost functions like path length or obstacle clearance [19]. The methodology was used to design a nonlinear compositional gradient in the Fe-Ni-Cr space that avoided the deleterious sigma and CrNi₂ phases at all relevant temperatures. In the study by Eliseeva et al. [20], this gradient was printed and, even with some compositional error during printing, shown to produce much less deleterious phases when compared with the linear gradient.

Though the previous methodology has been shown to successfully produce compositionally graded alloys without deleterious phases, it has not previously been able to optimize gradients for properties, despite the fact that the properties of a Functionally Graded Material are of critical importance to its design. While the properties of many polymer or metal-ceramic FGMs vary linearly or smoothly across the gradient region [1,2,5], compositionally graded alloys often have highly nonlinear property gradients. Properties also often change much more drastically and unpredictably with composition than with structure. Property gradients that are not obvious design objectives can also be important to the success of an FGM. For example, the discontinuities in stiffness and thermal expansion introduced by secondary phases are the very reason for their detriment to part integrity [15,16].

In this work, a new cost function is proposed that introduces property optimization into the previous FGM design methodology. It is argued that compositional paths with monotonic property profiles are preferred to all paths with non-monotonic profiles. It is further shown that monotonic property gradients can be made to take many diverse forms on a graded part by carefully controlling deposition rate. A metric for non-monotonicity from the literature [21] is then introduced into a cost function that is compatible with optimal path planners and finds the shortest feasible path with monotonic properties. A synthetic case study is then conducted to determine the effect of a cost function parameter on the balance between length and monotonicity. Additionally, a compositionally graded alloy is designed in the Fe-Co-Cr space to demonstrate the effectiveness of the methodology in finding a compositional path with a monotonic change in coefficient of thermal expansion (CTE). It is then shown how this compositional gradient could be deposited subject to a maximum deposition rate and property gradient.

In previous work [19,20], the authors have formulated the design of FGMs as a path planning problem consistent with motion planning frameworks [22,23]. The results of that formulation are summarized below.

The subset of undesirable compositions, like those with undesirable phases or with poor solidification properties, form the obstacle region, Z_obs ⊂ Z^d, that is to be avoided by the gradient path. In contrast, the complement of the obstacle region is the free space, $Zfree=Zd∖Zobs$⁠, which represents the viable composition space for gradient path planning.

In the previous work [19], a computational methodology was presented to solve the FGM design problem formulation and design optimal gradient paths. The methodology is reviewed below and visualized in Fig. 2.

calphad modeling [24–26] is used to predict the equilibrium phases present at a fixed thermodynamic state (i.e., composition and temperature). These predictions can be used to create an obstacle model or, in other words, a model of the locations of undesirable phase regions in composition space. While calphad could be sampled iteratively during the path planning process, most path planning algorithms requires millions of samples of the obstacle model to plan an optimal path. As such, the execution time of the obstacle model is of critical importance to the overall execution time of the path planning process. By using calphad samples to build a surrogate obstacle model, the computational expense of the obstacle model can be reduced by several orders of magnitude [19].

To build the surrogate obstacle model, a pseudo-random space-filling sampling (e.g., a Halton sequence sampling) of calphad is performed in the relevant composition space. The sampled compositions are then labeled as belonging to the obstacle region (e.g., those containing undesirable phases) or the free region. The labeled compositions are then given to a machine learning classifier which creates the surrogate obstacle model. The proper selection and training of a machine learning algorithm depends on the dimensionality of the composition space (i.e. the number of elements) as well as the complexity of the constraints. The purpose of the classifier is to minimize the computational cost needed to evaluate obstacle collision while maintaining an accurate representation of the obstacles.

In addition to the obstacle model, a property model might be necessary if the cost function is a function of properties. For example, in this work, a cost function is presented that prioritizes monotonic property profiles along a gradient path, but a property model is needed to predict these profiles. Suitable property models should be valid for the entire relevant composition space or at least the free region, Z_free. This condition could exclude many models for complex properties, but other, simpler models are valid for wide ranges of compositions like those for thermodynamic properties included in many calphad databases [26,27].

There is a wide selection of path planning algorithms that could solve the presented FGM design problem [28–31]. However, sampling-based planners were chosen given their adaptability to a wide variety of cost functions, scalability to high dimensions, and ease of implementation [22,30]. Specifically, Adiyatov and Varol’s [23] implementation of RRT* [22] is used in this work. Karaman and Frazzoli [22] developed RRT* to be an asymptotically optimal version of the Rapidly exploring Random Tree algorithm. Adiyatov and Varol [23] developed a fixed nodes implementation of RRT* (RRT*FN) that limits the maximum number of nodes in a tree and therefore limits the occupied computer memory.

In simple terms, RRT*FN iteratively samples the environment and creates a tree of paths in free space, Z_free. New samples are connected to a parent node in the tree so that the path formed minimizes the cost function, provided the connection would not collide with the obstacle region, Z_obs. Surrounding nodes are evaluated to determine if their paths would be cheaper with the new sample as their parent node and corresponding connections are rewired. Given enough random samples, RRT* is proven to find a feasible path if one exists (i.e., probabilistic completeness) [22]. Also, the optimal path in the tree approaches the globally optimal path as the number of samples increases (i.e., asymptotic optimality) [22].

Functionally Graded Materials are desired primarily for the gradients they produce in properties and therefore function [3–5]. As such, the consideration of properties should be a primary component of FGM design. Even properties that are not directly relevant to performance can be critical to manufacturability and part integrity. One such property is the CTE, which often varies dramatically with composition in alloys. In fact, large discontinuities in CTE are often present at phase boundaries in gradient materials. During the high-temperature additive manufacturing process, FGMs experience large thermal gradients. Consequently, these discontinuities in CTE often lead to cracking during manufacturing [15,16].

The compositional path that minimizes the stresses induced by property gradients is one that minimizes the property gradients themselves. Let y represent the position in physical space along a gradient part and p(z) represent some property of interest that is a function of material composition (e.g., CTE). If the local change in the property p along the physical part, dp/dy, is too large, then significant stress gradients can form. Therefore, the gradient path that minimizes property-induced stresses is that which minimizes dp/dy. The maximum property gradient dp/dy along a part, and consequently the maximum induced stress, is minimized when the property is graded linearly along the part and dp/dy is constant everywhere. However, it is unlikely that there exists a compositional gradient path σ(α) for which an arbitrary property varies linearly as a function of path index α and dp/dα is constant.

This condition mandates that properties vary monotonically along the compositional gradient path, meaning they either always decrease or always increase with path index, α. It also suggests that any property that is monotonic with path index can be mapped into a linear property gradient on the physical part by controlling material deposition rate, subject to the maximum and minimum achievable deposition rates. This process is depicted visually in Fig. 3.

Pathwise property profiles, p(α), that are monotonic with path index can not only be transformed to be linear along a part, they can also be mapped into any partwise property profile, p(y), that is monotonic with y. Furthermore, as illustrated in Fig. 4, pathwise property profiles that are monotonic with path index can even be transformed to vary non-monotonically along a part, given that the non-monotonic part profile is still bounded by the property values at z_init and z_goal. This capability means that a pathwise property profile, p(α), that is monotonic with path index can be mapped into nearly any properly bounded partwise profile, p(y). This characteristic is critical to the design of functionally graded parts. If a designer knows the desired maximum and minimum properties of a part, a material gradient with monotonic properties between those extrema can be used as a kind of “designer’s palette” of materials with which any spatial distribution of properties within the part can be created, subject to machine deposition rates.

By controlling deposition rates, gradient paths with monotonic property profiles can be mapped onto a part to achieve practically any spatial distribution of properties a designer might desire, subject to machine limitations. Therefore, to find a compositional gradient path that will provide a desired partwise property profile one need only find a path with properties that vary monotonically with path index. All feasible compositional gradient paths with monotonic property profiles offer this same design freedom and are therefore roughly equivalent in designer preference. Because non-monotonic profiles do not offer this freedom, all paths with monotonic profiles are preferred to those with non-monotonic profiles.

To introduce the ability to find gradient paths with monotonic properties into the current FGM design methodology, a cost function is needed that prioritizes such paths. Davydov and Zitikis [21] have developed a metric to quantify the degree of non-monotonicity possessed by a non-monotonic function. This metric, lack of monotonicity (LOM), relies on the calculation of two additional quantities: lack of increase (LOI) and lack of decrease (LOD).

To implement these metrics into the current methodology they must be formulated into a compatible cost function. RRT*FN [23] and many other tree- or map-based planners store segment costs independently and then sum the costs of each segment in a path to compute the total path cost. As such, a compatible cost function must be additive, satisfying the condition that c(σ₁|σ₂) = c(σ₁) + c(σ₂) for all σ₁, σ₂ ∈ Σ. This condition is incompatible with LOM because, as shown in Eq. (10), LOI and LOD must be evaluated for the total path before selecting which is the minimum quantity. For example, a path σ₁ with strictly increasing properties would have $LOMσ1=0$ and a path σ₂ with strictly decreasing properties would have $LOMσ2=0$⁠, but their concatenated path would have $LOM(σ1|σ2)≠0$⁠.

To remedy this issue, the change in properties between the initial and goal compositions can be evaluated and then either LOI or LOD could be used alone in place of LOM. For example, if Δp_total > 0 (i.e., p(z_init) < p(z_goal)), any monotonic profile between $p(zinit)$ and $p(zgoal)$ must be non-decreasing. In this case, LOI will measure the deviation from a monotonic property profile because all monotonic profiles must be non-decreasing (i.e., LOI = 0), while all non-monotonic profiles must increase somewhere (i.e., LOI > 0). Similarly, if Δp_total > 0 (i.e., p(z_init) > p(z_goal)), LOD will measure the deviation from all monotonic profiles, which are all non-increasing.

As reviewed previously, in order for a cost function to ensure asymptotic optimality with an optimal sampling-based planner like PRM*, RRT*, or RRT*FN, it must be strictly positive, monotonic, and bounded [22,23]. Both LOI and LOD are bounded by $∫ab|g0′|dλ$ [21] and monotonic for concatenated paths, but neither are strictly positive. This is because $LOIσ(p)=0$ and $LODσ(p)=0$ for non-decreasing and non-increasing property profiles, respectively. As such, there is no way for LOI or LOD to discriminate between monotonic property profiles.

However, by including path length (i.e., the Euclidean distance traversed by the path in composition space) into the cost function with LOI or LOD, shorter monotonic paths will become preferred to longer ones. Also, the cost function itself becomes strictly positive as path length is strictly positive and both LOI and LOD are non-negative. Designers can diminish the effect of path length in the cost function by introducing a weighting parameter, w, with very low magnitude. Such a parameter could effectively ensure path length is an active objective only when LOI or LOD is zero. A study on the magnitude of this parameter is presented later in this work.

The cost function shown in Eq. (13) was designed to seek the shortest gradient path with a monotonic property profile. A simple synthetic case study was created to test the ability of the cost function to find such paths. The case study was also used to examine the effect of the parameter w on the cost function’s tendency to balance length and lack of monotonicity.

An artificial property model was generated for a two-dimensional input space, x₁, x₂ ∈ [0, 1]. The property model was created to have non-monotonic regions, but still enable monotonic paths from z_init = (1, 1) to z_goal = (0, 0). This model, p(x₁, x₂), was created from an initial planar surface rotated along the line x₁ + x₂ = 1 to be monotonic with both x₁ and x₂. A non-monotonic region was created on the planar surface by modeling a semi-ellipsoid on the surface of the plane that was oriented normal to the plane with a major axis parallel to the line x₁ + x₂ = 1. The semi-ellipsoid was sized to be small enough to allow for monotonic paths from z_init = (1, 1) to z_goal = (0, 0), but large enough to make the straight-line, shortest length path non-monotonic. The final model has values ranging from approximately 0.3 to 0.7 for x₁, x₂ ∈ [0, 1]. Figure 5 illustrates the surface of the synthetic property model.

To test the cost function presented in Eq. (13), the path planning algorithm, RRT*FN, was used to plan several paths in x₁ and x₂ from z_init = (1, 1) to z_goal = (0, 0). For simplicity, no obstacle region was considered. Because the synthetic property model has a lower value at z_goal than z_init, LOD was used to measure non-monotonicity. Each of the seven runs of the path planning algorithm used a different value of the parameter w in the cost function, increasing in magnitude from 10⁻⁶ to 1. As such, each run had a different weighting between lack of decrease and path length. Each run of RRT*FN generated 5000 random samples with the same random seed, meaning each run used the same points to generate a tree and find an optimal path.

The paths planned by each run of RRT*FN are shown in Fig. 6, and their respective pathwise property profiles are shown in Fig. 7. The values of each term in the cost function for the optimal paths in each run are listed in Table 1. When the weighting parameter is very small (w ≤ 10⁻³), the cost function seems to perform as intended and prioritize monotonicity before path length. In fact, the paths produced when $w=10−6,10−5,10−4,and10−3$ are exactly the same and resemble the desired path: the shortest monotonic path between z_init and z_goal. As w increases to 10⁻² and 10⁻¹, the paths become slightly shorter by encroaching into the ellipsoidal region, but also produce non-monotonic property profiles (⁠$LODσ(p)>0$⁠). When w is further increased to 1, the length term in the cost function completely dominates and the path produced is the straight-line, shortest length path between z_init and z_goal.

In practice, w should be set as low as possible to make minimizing length a secondary objective to promoting monotonicity. However, these results indicate that a value of w ≤ 10⁻³ should produce the desired behavior for similarly scaled problems. All composition spaces are bounded by 0 and 1, so most FGM design problems should be of similar scale if the property model is scaled to have bounds near 0 and 1.

A second case study was conducted to examine the performance of the proposed cost function in a realistic FGM design problem. Specifically, a compositionally graded alloy was designed with a monotonic gradient in CTE. As discussed earlier, steep gradients in CTE can often lead to detrimental stress gradients during the additive manufacturing process [15,16] and are therefore of critical importance to all compositionally graded alloys.

The iron-cobalt-chromium (Fe-Co-Cr) ternary system was chosen for this case study because it is easily visualized in two dimensions and has significant relevance to many engineering materials. It is important to note that the methodology can be applied to a higher dimensional system with more elements, but a simpler system was chosen to visualize how optimal paths navigate the space. The end points of the designed gradient path, z_init and z_goal, were chosen to be Fe95Co5 [at%] and Fe10Co60Cr30 [at%], respectively. The initial point, Fe95Co5 [at%], is representative of many steels including high-speed steels used in cutting tools for their high temperature resistance and hardness [32]. The goal point, Fe10Co60Cr30 [at%], represents cobalt-chrome alloys that exhibit exceptional corrosion, wear and thermal resistance as well as high specific strength and are often used in biomedical applications [33]. A gradient between these two materials could conceivably be employed in a surgical device or a device with significant thermal requirements like a drill or turbine.

To apply the FGM design methodology to this system, phase regions were first modeled with calphad software, specifically Thermo-Calc’s TCHEA2 database [26,27]. The calphad model was sampled 1275 times in a regular grid throughout the Fe-Co-Cr composition space. Phase equilibria were calculated at a temperature of 1000 K, which was chosen to approximate the temperature of the manufacturing process. Sigma (σ) phase was chosen as an undesirable phase due to its detrimental mechanical properties including high brittleness [34,35]. A surrogate obstacle model was created by labeling the samples with greater than one percent mole fraction of sigma phase and then training a k-nearest neighbors (k = 3) classifier to represent the sigma phase region.

Thermo-Calc’s TCHEA2 can also predict CTE for any composition in its database [27]. This capability was leveraged to gather predictions of CTE at 1000 K at the same 1275 compositions that were sampled for phase information. The CTE data was made to lie within 1 × 10⁻⁶ K⁻¹ and 1 × 10⁻⁴K⁻¹ by rounding outlying data to the bounds and was then scaled to range from 0 to 1. An interpolant model was then created from the data that linearly interpolates between the data to predict CTE at any composition in the Fe-Co-Cr composition space.

Given a model of the obstacle region and a model of the relevant property, the path planning algorithm, RRT*FN [23], was used to plan paths for two different cost functions. The first cost function was simply path length and was planned to assess the properties of a path planned without any consideration of optimal properties. The second cost function was that proposed in Eq. (13), which seeks to find the shortest path with monotonic properties. Because the CTE at z_init is greater than that at z_goal, $LODσ(p)$ was chosen to measure non-monotonicity. As in the synthetic case study, the random seed was fixed for both cases so the same 5000 nodes were used to construct the tree and find the optimal path.

The paths found to optimize each cost function are displayed in Fig. 8 and the value of CTE along each path is shown in Fig. 9. The path planned for the first cost function, c(σ) = l, resembles the straight-line path between z_init and z_goal. While this is the shortest feasible path, it experiences a large increase in CTE as it approaches z_goal, as seen in Fig. 9. This makes the property profile significantly non-monotonic. The second cost function however, succeeds in producing the shortest path with a monotonic CTE profile.

While the path planned with the proposed cost function is monotonic, Fig. 9 does show a sharp decrease in CTE. However, the shortest length path also experiences a similar drop in CTE. In fact, by examining the CTE contours in Fig. 8, it appears that such a drop is unavoidable as every feasible path to z_goal would experience a steep decrease in CTE. Fortunately, these steep regions can be mitigated by decreasing the deposition rate in these regions. Although the first path is shorter, it has more steep changes in CTE and might actually need more layers to print due to the need to decrease deposition rates in steep regions. In addition to potentially requiring less layers to print, a part made with the monotonic compositional gradient could be made to have nearly any properly bounded CTE profile, subject to machine limitations on deposition rate, as shown in Fig. 4.

To demonstrate how one could determine deposition rates, the deposition of the monotonic path was planned for a hypothetical part. While the path could be deposited so that the partwise property, p(y), was linear, it would require extreme differences in the deposition rates because of the differences in dp/dα between the steep drop in CTE and the comparatively flat regions at the beginning and end of the path, as seen in Fig. 9. This would lead to a relatively enormous region of the part dedicated to a small portion of the compositional gradient.

Instead, reasonable constraints for the partwise property gradient and deposition rate were used to determine the deposition plan. The position on the part, y, was measured in build layers for simplicity. The maximum property gradient on the part, $|dp/dy|max$⁠, was chosen to be 4 × 10⁻⁷K⁻¹ per layer. The achievable change in composition per layer was estimated to be 1 atomic percent [at%]. This number was then divided by the total length of the path to estimate a maximum deposition rate $|dα/dy|max$⁠. These two constraints were then used with Eq. (6) to determine how many layers should be used to print each segment of the path.

Figure 10 displays the resulting gradients of path index, α, and CTE along the planned gradient part. Note that there are three distinct regions where a different constraint was active. In the first part of the gradient, where dp/dα is small, the deposition rate is at its maximum. This is where the gradient is deposited as fast as possible because the associated property gradient dp/dy is insignificant. Once the steep decline in CTE is reached, the deposition rate is slowed for about 100 layers so that the property gradient is at its maximum acceptable value, $|dp/dy|max$⁠. After the steep change in CTE has been navigated, the deposition rate is once again maximized. The compositions of the corresponding part are shown in Fig. 11. From layers 100 to 200, where deposition rate is slowest, the compositions barely change. This is to diminish the effects of the steep property gradient, dp/dα, seen in Fig. 9.

In prior work, a computational design methodology was developed to plan compositionally graded alloys that avoid undesirable phase regions [19]. This methodology is extended by introducing a new cost function with a metric from literature [21] that measures non-monotonicity.With the improved methodology, designers can now find the shortest feasible gradient path with a monotonic property profile.

By controlling material deposition rate, compositional paths with monotonic property profiles can be made into gradient parts with nearly any properly bounded spatial distribution of properties. Consequently, a material gradient with a monotonic change in properties can be used as a kind of “designer’s palette” of materials, wherein a one-to-one correspondence of property to material can be exploited to design properties locally. Because of this design freedom, designing a single gradient with monotonic properties avoids the need to design a new material gradient for every property gradient a designer might desire.

When using a gradient with monotonic properties, controlling the spatial distribution of properties is simply a matter of controlling material deposition. In this work, a straightforward, one-dimensional deposition strategy was employed to satisfy constraints on maximum deposition rate and maximum property gradient. However, more complex deposition strategies can be easily designed to achieve desired property distributions by inverting the monotonic, and therefore one-to-one, material-property relationship.

Techniques like Multi-Material Topology Optimization (MMTO) [36–38], can be used to optimize the material distribution, as well as the topology, of a functionally graded part. To the authors knowledge, these methods have not been applied to compositionally graded alloys, but have been applied to linear gradients in composite or polymer materials with simpler property profiles [36,38]. Because the property gradients produced by the proposed cost function are monotonic, these methods could likely be adapted to optimize material distribution in compositionally graded metal parts by using path index as a design variable instead of volume fraction. Such an integration with MMTO techniques would enable a full material-to-part design process for functionally graded metal parts, but is left to future work.

FGMs are often used for their ability to achieve incompatible performance objectives that a single monolithic material could not. Although the proposed cost function provides a means to optimize gradient paths for a single property, it has not yet been shown to optimize multiple properties at once. It is not clear, for example, how best to optimize two inversely correlated property gradients in a single FGM as both properties will be tightly coupled and some preference between each property might need to be established. The ability to design gradients for multiple properties would be an important contribution, but the development of such a framework is also left to future work.

Additionally, there are multiple sources of uncertainty in the proposed methodology including the calphad thermodynamic models, the property models, and the surrogate modeling techniques. Evolving the methodology to plan robust gradient paths that consider these uncertainties is yet another avenue for future research.

T. K. acknowledges the support of NSF through the project NRT-DESE: Data-Enabled Discovery and Design of Energy Materials (D3EM), NSF-DGSE- 1545403.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

c =

cost function that relates a path σ to a measure of preference

d =

number of relevant elements in the functionally graded material

l =

length of a path in state space

p =

material property that is dependent on composition and, therefore, position on a gradient part

r =

the ratio of the length of path σ₁ to the length of the concatenation of paths σ₁ and σ₂

w =

weighting parameter to control the effect of path length in the proposed cost function

y =

position on a gradient part

z =

point in the state space of the path planning problem, Z^d

x_i =

composition of element i

z_goal =

goal point in path (final composition)

z_init =

initial point in path (initial composition)

Z_free =

obstacle-free region (region in the composition space with allowable characteristics)

Z_obs =

obstacle region (region in the composition space with unacceptable characteristics)

Z^d =

state space of the path planning problem (composition space of a functionally graded material)

α =

index of a point in the path (the path length up to the indexed position normalized by the total path length)

Δy_i→i+1,min =

the minimum length in physical space (i.e., number of layers) needed to print an arbitrary segment of a gradient, subject to constraints

σ =

continuous path in the composition space

σ_best =

feasible path that minimizes the cost function c