## Introduction

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 [35]. 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 [35].

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 [69]. 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 [1116]. 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 CrNi2 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.

## Methods

### Functionally Graded Material Design as a Path Planning Problem.

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.

Let Zd represent the relevant composition space of the gradient, where d is the dimensionality of the space and is equal to the number of relevant elements. Let z represent a point in that space: the composition of every relevant element i at a single material point, as seen in Eq. (1).
$z={x1,…,xd:∑i=1dxi=1andxi≥0∀i}$
(1)

The subset of undesirable compositions, like those with undesirable phases or with poor solidification properties, form the obstacle region, ZobsZd, 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.

The goal of the FGM design problem is to find the optimal gradient or path in this composition space between two predefined compositions, zinit and zgoal. Let the continuous function σ : [0, 1] → z represent such a path. Let α ∈ [0, 1] represent a path index that scales with distance traveled along a path (i.e., path length), l, in state space (i.e., composition space). Equation (2) represents the scaling of α where α is defined as the path length up to a position σ(α) along the path normalized by the total length of the path.
$α:=l(σ(0),σ(α))l(σ(0),σ(1))∀α∈[0,1]$
(2)
A compositional gradient path is collision-free if it is contained entirely within free space: σ(α) ∈ Zfree$∀$α ∈ [0, 1]. A gradient path is feasible if it is collision free, σ(0) = zinit, and σ(1) = zgoal. Figure 1 illustrates a feasible path in a simplified example. The optimal gradient path, σbest, is that which is feasible and also minimizes the cost function c : σR≥0. The cost function (e.g., path length) is defined by the designer, but must be strictly positive. To be strictly positive, the cost function must equal zero if and only if $σ(α)=σ(0),∀α∈[0,1]$. The FGM design problem is summarized in the problem formulation below.
$Findσbest=argminσc(σ)subjecttoσ(α)∈Zfree∀α∈[0,1],σ(0)=zinit,σ(1)=zgoal$

### Computational Functionally Graded Material Design Methodology.

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 [2426] 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, Zfree. 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 [2831]. 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, Zfree. 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, Zobs. 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].

As discussed by Karaman and Frazzoli [22], for any optimal sampling-based path planning algorithm (e.g., PRM*, RRT*, RRT*FN) to be asymptotically optimal, cost functions must be strictly positive, monotonic, and bounded. To be strictly positive, the cost of a path can only be zero if the path never moves from its starting position (i.e., c(σ) = 0 if and only if $σ(α)=σ(0),∀α∈[0,1]$). A monotonic cost function satisfies the condition that c(σ1) ≤ c(σ1|σ2) for all σ1, σ2 ∈ Σ, where σ1|σ2 represents the concatenation of paths. The concatenation of paths is an operation performed on paths σ1, σ2 ∈ Σ where σ1(1) = σ2(0). The concatenated path is defined by Eq. (3), where $r=lσ1/(lσ1+lσ2)$ is the ratio between the length of path σ1 and the combined lengths of paths σ1 and σ2.
$(σ1|σ2)(α):={σ1(αr)∀α∈[0,r]σ2(α−r1−r)∀α∈[r,1]$
(3)

### Property Monotonicity as a Design Objective.

Functionally Graded Materials are desired primarily for the gradients they produce in properties and therefore function [35]. 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.

Fortunately, the rate at which the compositional gradient is deposited on the physical part is an accessible design parameter that can be varied during the build. Consider dα/dy to be the local change in path index per unit length along a gradient part (i.e., per build layer). Because path index scales proportionally with length in composition space (Eq. (2)) and each path index directly corresponds to a specific material composition, this rate is analogous to deposition rate: the rate that material composition is changing per unit length along a gradient part. This rate can be approximated by measuring the composition of two adjacent layers in a gradient part and dividing the difference in the associated path indices by the layer thickness. Controlling this deposition rate could allow designers to vary properties linearly along the dimensions of the part even if the property does not vary linearly with composition or path index. For example, due to the relationship shown in Eq. (4), one could achieve a constant dp/dy by adjusting dα/dy in inverse proportion with dp/dα as α increases from 0 to 1.
$dpdy=dpdα(dαdy)$
(4)
However, this is only possible if the condition shown in Eq. (5) is true, where Δptotal = p(zgoal) − p(zinit).
$sgn(dpdα)=sgn(Δptotal)∀α∈[0,1]$
(5)

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 zinit and zgoal. 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.

To determine the deposition plan needed to achieve a desired partwise property profile, one need only manipulate the relationship shown in Eq. (4). In practice, the deposition of a gradient path might be designed such that practical constraints on the maximum property gradient $|dp/dy|max$ and maximum deposition rate $|dα/dy|max$ are satisfied. Equation (6) provides a method for determining the deposition plan for a gradient path subject to these constraints. Given the function p(α) that describes the change in the relevant property as a function of path index, one can divide the path in to segments where p(α) is approximately linear. For an arbitrary segment between αi and αi+1, Eq. (6) describes the minimum length of the gradient part (i.e., the minimum number of layers) needed to print that segment of the gradient path, Δyii+1,min.
$Δyi→i+1,min=max{|αi+1−αi||dαdy|max−1,|p(αi+1)−p(αi)||dpdy|max−1}$
(6)
Equation (6) can be calculated for each segment of the path to construct a piecewise linear deposition plan y(α), as shown in Eq. (7), for which every position or layer on a gradient part has an associated path index and therefore material composition.
$y(αi+1)=y(αi)+Δyi→i+1,min$
(7)

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.

### Lack of Monotonicity.

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 calculate the lack of increase of a general continuous function g0 over the interval [a, b] one must integrate the negative part of the first derivative of g0 over the Lebesgue measure, λ, as seen in Eq. (8) where $(g0′)−=max{−g0′,0}$. Note that if g0 is monotonically increasing over the interval [a, b], then LOI[a,b](g0) = 0.
$LOI[a,b](g0)=∫ab(g0′)−dλ$
(8)
Lack of decrease is calculated similarly by integrating the positive part of the first derivative of g0, as seen in Eq. (9) where $(g0′)+=max{g0′,0}$.
$LOD[a,b](g0)=∫ab(g0′)+dλ$
(9)
Finally, the lack of monotonicity of g0 over the interval [a, b] is shown in Eq. (10).
$LOM[a,b](g0)=2min{LOI[a,b](g0),LOD[a,b](g0)}$
(10)
These metrics can be used to assess the lack of monotonicity in a pathwise property profile and consequently drive the FGM design methodology toward gradient paths with monotonic properties. The LOI and LOD of the property profile, p(α), for a given path, σ, can be calculated from Eqs. (11) and (12), respectively.
$LOIσ(p)=∫01(dpdα)−dα$
(11)
$LODσ(p)=∫01(dpdα)+dα$
(12)

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(σ1|σ2) = c(σ1) + c(σ2) for all σ1, σ2 ∈ Σ. 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 σ1 with strictly increasing properties would have $LOMσ1=0$ and a path σ2 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 Δptotal > 0 (i.e., p(zinit) < p(zgoal)), 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 Δptotal > 0 (i.e., p(zinit) > p(zgoal)), 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.

Equation (13) presents a cost function that satisfies all the aforementioned criteria. As such, it can be easily implemented into the current methodology to obtain gradient paths with monotonic property profiles.
$c(σ)={LOIσ(p)+wlifΔptotal>0LODσ(p)+wlifΔptotal<0$
(13)

## Case Study Results and Discussion

### Synthetic Case Study.

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, x1, x2 ∈ [0, 1]. The property model was created to have non-monotonic regions, but still enable monotonic paths from zinit = (1, 1) to zgoal = (0, 0). This model, p(x1, x2), was created from an initial planar surface rotated along the line x1 + x2 = 1 to be monotonic with both x1 and x2. 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 x1 + x2 = 1. The semi-ellipsoid was sized to be small enough to allow for monotonic paths from zinit = (1, 1) to zgoal = (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 x1, x2 ∈ [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 x1 and x2 from zinit = (1, 1) to zgoal = (0, 0). For simplicity, no obstacle region was considered. Because the synthetic property model has a lower value at zgoal than zinit, 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−6 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−3), 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 zinit and zgoal. As w increases to 10−2 and 10−1, 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 zinit and zgoal.

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−3 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.

### Materials Case Study.

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, zinit and zgoal, 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−6 K−1 and 1 × 10−4K−1 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 zinit is greater than that at zgoal, $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 zinit and zgoal. While this is the shortest feasible path, it experiences a large increase in CTE as it approaches zgoal, 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 zgoal 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−7K−1 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.

## Conclusions and Future Work

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) [3638], 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.

## Acknowledgment

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

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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

## Nomenclature

• 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 σ1 to the length of the concatenation of paths σ1 and σ2

•
• w =

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

•
• y =

•
• z =

point in the state space of the path planning problem, Zd

•
• xi =

composition of element i

•
• zgoal =

goal point in path (final composition)

•
• zinit =

initial point in path (initial composition)

•
• Zfree =

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

•
• Zobs =

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

•
• Zd =

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)

•
• Δyii+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

## References

References
1.
Koizumi
,
M.
, and
Niino
,
M.
,
1995
, “
Overview of FGM Research in Japan
,”
MRS Bull.
,
20
(
1
), pp.
19
21
. 10.1557/S0883769400048867
2.
Markworth
,
A. J.
,
Ramesh
,
K. S.
, and
Parks
,
W. P.
,
1995
, “
Modelling Studies Applied to Functionally Graded Materials
,”
J. Mater. Sci.
,
30
(
9
), pp.
2183
2193
. 10.1007/BF01184560
3.
Oxman
,
N.
,
2011
, “
Variable Property Rapid Prototyping
,”
Virtual Phys. Prototyping
,
6
(
1
), pp.
3
31
. 10.1080/17452759.2011.558588
4.
Vaezi
,
M.
,
Chianrabutra
,
S.
,
Mellor
,
B.
, and
Yang
,
S.
,
2013
, “
Multiple Material Additive Manufacturing—Part 1: A Review
,”
Virtual Phys. Prototyping
,
8
(
1
), pp.
19
50
. 10.1080/17452759.2013.778175
5.
Loh
,
G. H.
,
Pei
,
E.
,
Harrison
,
D.
, and
Monzón
,
M. D.
,
2018
, “
,”
,
23
, pp.
34
44
6.
Hofmann
,
D. C.
,
Kolodziejska
,
J.
,
Roberts
,
S.
,
Otis
,
R.
,
Dillon
,
R. P.
,
Suh
,
J. -O.
,
Liu
,
Z. -K.
, and
Borgonia
,
J. -P.
,
2014
, “
,”
J. Mater. Res.
,
29
(
17
), pp.
1899
1910
. 10.1557/jmr.2014.208
7.
Tammas-Williams
,
S.
, and
Todd
,
I.
,
2017
, “
Design for Additive Manufacturing With Site-specific Properties in Metals and Alloys
,”
Scr. Mater.
,
135
, pp.
105
110
. 10.1016/j.scriptamat.2016.10.030
8.
Yan
,
L.
,
Chen
,
Y.
, and
Liou
,
F.
,
2020
, “
,”
,
31
, p.
100901
9.
Reichardt
,
A.
,
Shapiro
,
A. A.
,
Otis
,
R.
,
Dillon
,
R. P.
,
Borgonia
,
J. P.
,
McEnerney
,
B. W.
,
Hosemann
,
P.
, and
Beese
,
A. M.
,
2020
, “
,”
Int. Mater. Rev.
, pp.
1
29
. 10.1080/09506608.2019.1709354
10.
Schwendner
,
K. I.
,
Banerjee
,
R.
,
Collins
,
P. C.
,
Brice
,
C. A.
, and
Fraser
,
H. L.
,
2001
, “
Direct Laser Deposition of Alloys From Elemental Powder Blends
,”
Scr. Mater.
,
45
(
10
), pp.
1123
1129
. 10.1016/S1359-6462(01)01107-1
11.
Hofmann
,
D. C.
,
Roberts
,
S.
,
Otis
,
R.
,
Kolodziejska
,
J.
,
Dillon
,
R. P.
,
Suh
,
J.-o.
,
Shapiro
,
A. A.
,
Liu
,
Z. -K.
, and
Borgonia
,
J. -P.
,
2014
, “
,”
Sci. Rep.
,
4
, p.
5357
. 10.1038/srep05357
12.
Carroll
,
B. E.
,
Otis
,
R. A.
,
Borgonia
,
J. P.
,
Suh
,
J.-o.
,
Dillon
,
R. P.
,
Shapiro
,
A. A.
,
Hofmann
,
D. C.
,
Liu
,
Z. -K.
, and
Beese
,
A. M.
,
2016
, “
Functionally Graded Material of 304L Stainless Steel and Inconel 625 Fabricated by Directed Energy Deposition: Characterization and Thermodynamic Modeling
,”
Acta. Mater.
,
108
, pp.
46
54
. 10.1016/j.actamat.2016.02.019
13.
Chen
,
B.
,
Su
,
Y.
,
Xie
,
Z.
,
Tan
,
C.
, and
Feng
,
J.
,
2020
, “
Development and Characterization of 316L/Inconel625 Functionally Graded Material Fabricated by Laser Direct Metal Deposition
,”
Optics Laser Technol.
,
123
, p.
105916
. 10.1016/j.optlastec.2019.105916
14.
Reichardt
,
A.
,
Dillon
,
R. P.
,
Borgonia
,
J. P.
,
Shapiro
,
A. A.
,
McEnerney
,
B. W.
,
Momose
,
T.
, and
Hosemann
,
P.
,
2016
, “
Development and Characterization of Ti-6Al-4V to 304L Stainless Steel Gradient Components Fabricated With Laser Deposition Additive Manufacturing
,”
Mater. Des.
,
104
, pp.
404
413
. 10.1016/j.matdes.2016.05.016
15.
Bobbio
,
L. D.
,
Otis
,
R. A.
,
Borgonia
,
J. P.
,
Dillon
,
R. P.
,
Shapiro
,
A. A.
,
Liu
,
Z. -K.
, and
Beese
,
A. M.
,
2017
, “
Additive Manufacturing of a Functionally Graded Material From Ti-6Al-4V to Invar: Experimental Characterization and Thermodynamic Calculations
,”
Acta. Mater.
,
127
, pp.
133
142
. 10.1016/j.actamat.2016.12.070
16.
Meng
,
W.
,
Xiaohui
,
Y.
,
Zhang
,
W.
,
Junfei
,
F.
,
Lijie
,
G.
,
Qunshuang
,
M.
, and
Bing
,
C.
,
2020
, “
Additive Manufacturing of a Functionally Graded Material From Inconel625 to Ti6Al4V by Laser Synchronous Preheating
,”
J. Mater. Process. Technol.
,
275
, p.
116368
. 10.1016/j.jmatprotec.2019.116368
17.
Moustafa
,
A. R.
,
Durga
,
A.
,
Lindwall
,
G.
, and
Cordero
,
Z. C.
,
2020
, “
Scheil Ternary Projection (STeP) Diagrams for Designing Additively Manufactured Functionally Graded Metals
,”
,
32
, p.
101008
18.
Bocklund
,
B.
,
Bobbio
,
L. D.
,
Otis
,
R. A.
,
Beese
,
A. M.
, and
Liu
,
Z.-K.
,
2020
, “
,”
Materialia
,
11
, p.
100689
. 10.1016/j.mtla.2020.100689
19.
Kirk
,
T.
,
Galvan
,
E.
,
Malak
,
R.
, and
Arroyave
,
R.
,
2018
, “
,”
ASME J. Mech. Des.
,
140
(
11
), p.
111410
. 10.1115/1.4040816
20.
Eliseeva
,
O. V.
,
Kirk
,
T.
,
Samimi
,
P.
,
Malak
,
R.
,
Arróyave
,
R.
,
Elwany
,
A.
, and
Karaman
,
I.
,
2019
, “
Functionally Graded Materials Through Robotics-inspired Path Planning
,”
Mater. Des.
,
182
, p.
107975
. 10.1016/j.matdes.2019.107975
21.
Davydov
,
Y.
, and
Zitikis
,
R.
,
2017
, “
Quantifying Non-monotonicity of Functions and the Lack of Positivity in Signed Measures
,”
Modern Stochastics: Theory and Applications
,
4
(
3
), pp.
219
231
.
arXiv: 1705.02742
. 10.15559/17-VMSTA84
22.
Karaman
,
S.
, and
Frazzoli
,
E.
,
2011
, “
Sampling-Based Algorithms for Optimal Motion Planning
,”
Int. J. Rob. Res.
,
30
(
7
), pp.
846
894
. 10.1177/0278364911406761
23.
,
O.
, and
Varol
,
H.
,
2013
, “
Rapidly-Exploring Random Tree Based Memory Efficient Motion Planning
,”
2013 IEEE International Conference on Mechatronics and Automation (ICMA)
,
Takamatsu, Japan
, IEEE, pp.
354
359
.
24.
Kaufman
,
L.
, and
Bernstein
,
H.
,
1970
,
Computer Calculation of Phase Diagrams. With special Reference to Refractory Metals
,
,
New York
.
25.
Kaufman
,
L.
, and
Agren
,
J.
,
2014
, “
CALPHAD, First and Second Generation–Birth of the Materials Genome
,”
Scr. Mater.
,
70
, pp.
3
6
. 10.1016/j.scriptamat.2012.12.003
26.
,
J.-O.
,
Helander
,
T.
,
Höglund
,
L.
,
Shi
,
P.
, and
Sundman
,
B.
,
2002
, “
Thermo-Calc & DICTRA, Computational Tools for Materials Science
,”
,
26
(
2
), pp.
273
312
. 10.1016/S0364-5916(02)00037-8
27.
Mao
,
H.
,
Chen
,
H.-L.
, and
Chen
,
Q.
,
2017
, “
TCHEA1: A Thermodynamic Database Not Limited for “High Entropy” Alloys
,”
J. Phase. Equilib. Diffus.
,
38
(
4
), pp.
353
368
. 10.1007/s11669-017-0570-7
28.
Hwang
,
Y. K.
, and
Ahuja
,
N.
,
1992
, “
Gross Motion Planning—A Survey
,”
ACM Comput. Surveys (CSUR)
,
24
(
3
), pp.
219
291
. 10.1145/136035.136037
29.
Raja
,
P.
, and
Pugazhenthi
,
S.
,
2012
, “
Optimal Path Planning of Mobile Robots: A Review
,”
Int. J. Phys. Sci.
,
7
(
9
), pp.
1314
1320
. 10.5897/IJPS11.1745
30.
Elbanhawi
,
M.
, and
Simic
,
M.
,
2014
, “
Sampling-based Robot Motion Planning: A Review
,”
IEEE Access
,
2
, pp.
56
77
. 10.1109/ACCESS.2014.2302442
31.
Mac
,
T. T.
,
Copot
,
C.
,
Tran
,
D. T.
, and
De Keyser
,
R.
,
2016
, “
Heuristic Approaches in Robot Path Planning: A Survey
,”
Rob. Auto. Syst.
,
86
, pp.
13
28
. 10.1016/j.robot.2016.08.001
32.
Gulyaev
,
A. P.
, and
Kupalova
,
I. K.
,
1970
, “
Effect of Cobalt on the Structure and Properties of High-Speed Steels
,”
Met. Sci. Heat. Treat.
,
12
(
8
), pp.
666
671
. 10.1007/BF00654791
33.
Kereiakes
,
D. J.
,
Cox
,
D. A.
,
Hermiller
,
J. B.
,
Midei
,
M. G.
,
Bachinsky
,
W. B.
,
Nukta
,
E. D.
,
Leon
,
M. B.
,
Fink
,
S.
,
Marin
,
L.
, and
Lansky
,
A. J.
,
2003
, “
Usefulness of a Cobalt Chromium Coronary Stent Alloy
,”
Am. J. Cardiology
,
92
(
4
), pp.
463
466
. 10.1016/S0002-9149(03)00669-6
34.
Hou
,
J.
,
Guo
,
J.
,
Zhou
,
L.
, and
Ye
,
H.
,
2006
, “
Sigma Phase Formation and Its Effect on Mechanical Properties in the Corrosion-Resistant Superalloy K44
,”
Z. Metallkunde
,
97
(
2
), pp.
174
181
. 10.3139/146.101222
35.
Hsieh
,
C.-C.
, and
Wu
,
W.
,
2012
, “
Overview of Intermetallic Sigma Phase Precipitation in Stainless Steels
,”
ISRN Metallurgy
,
2012
. 10.5402/2012/732471
36.
Xia
,
Q.
, and
Wang
,
M. Y.
,
2008
, “
Simultaneous Optimization of the Material Properties and the Topology of Functionally Graded Structures
,”
Comput.-Aided Design
,
40
(
6
), pp.
660
675
37.
Hiller
,
J. D.
, and
Lipson
,
H.
,
2009
, “
Multi Material Topological Optimization of Structures and Mechanisms
,”
Proceedings of the 11th Annual Conference on Genetic and Evolutionary Computation, GECCO ’09
,
ACM
,
, pp.
1521
1528
.
38.
Garland
,
A.
, and
,
G.
,
2015
, “
Design and Manufacturing Functionally Gradient Material Objects With An Off the Shelf Three-Dimensional Printer: Challenges and Solutions
,”
ASME J. Mech. Des.
,
137
(
11
), p.
111407
. 10.1115/1.4031097