A Dynamic Programming Setting for Functionally Graded Thick-Walled Cylinders

Material property variation in non-homogeneous internally pressurized thick-walled cylinders is investigated within the context of dynamic programming theory. The material is assumed to be linear, elastic, isotropic, and functionally graded in the radial direction. Based on the plane stress hypothesis, a state space formulation is given and the optimal control problem is stated and solved by means of Pontryagin’s Principle for different objective functionals. Optimal Young’s modulus distribution is found to be piecewise linear along the radial domain. A brief digression on the possible existence of switching points is addressed. Finally, a numerical example is performed within a special class of derived optimal solutions, showing promising results in terms of equivalent stress reduction with respect to the most used variations in literature.


Introduction
Functionally Graded (FG) materials are a kind of composite material whose physical and mechanical properties vary spatially along specific directions over the entire domain. They are present in many engineering applications such as, for example, space planes, nuclear fusion reactors, energy conversion systems, and thermo generators [1]. Various problems of FG materials have attracted considerable attention in recent years and their increasing use needs comprehensive mechanical analyses. A significant amount of work has been devoted to understand the thermal and mechanical behavior of FG axisymmetric bodies. In particular, the study of the stress and strain fields in FG cylinders [2,3], spheres [4,5], pressure vessels [6,7], circular plates [8,9], and rotating disks [10][11][12] has been widely investigated. In all the above-mentioned works, a generic property (e.g., Young's modulus, mass density, or a thermal expansion coefficient) is assumed to vary along the radial coordinate according to some pre-defined functions. Among all possible variations, power-law distributions are widely used [13]. A few works consider different property distribution laws. For instance, in [14], both exponential and parabolic variations of the mechanical properties are considered. Moreover, in [15], material properties are assumed to be exponential functions and exact solutions are obtained by the use of hypergeometric functions.
The optimum response of material properties to an actual environment is the main requirement in the design of FG materials; nevertheless, very few articles deal with the problem of optimal material distribution because, often, the optimization process heavily relies on subsequent Finite Element (FE) simulations. However, when the governing differential equations admit a closed form solution, which is the case of the power-law distributions, a different approach can be followed (see [16,17], for instance). In these cases, in fact, the decision variables in the optimization problem appear as parameters in the expression of the solution of the elastic problem. Consequently, the optimization process turns out to be equivalent to the minimization of a function, i.e., a static optimization problem.
Despite the fact that, when mechanical properties vary according to a power-law distribution, an analytic solution of the governing differential equations can be obtained, a more realistic (and, in some sense, intrinsic) optimization problem should consider all possible classes of property variations. Correspondingly, the solution of the optimization problem will be a function (to be chosen in a wider set), not just a set of parameters to be used as tuning values for a prefixed class of functions. To the extent of our knowledge, no previous works deal with the problem of finding optimal property distributions in FG axisymmetric mechanical components. Herein, a property variation is referred to as optimal if it maximizes or minimizes a (pre-defined) mechanical performance. For instance, a designer may be interested in the search for mechanical or physical properties variation for a vessel of minimum compliance, a cylinder of minimum hoop stress at a specific radius or a disk of uniform strength or minimum mass, complying with other constraints which can be economic or related to the manufacturing process. These problems can be formulated either in the context of optimal control theory (i.e., in a dynamic programming setting) or by means of calculus of variations. However, the former is more suitable when the optimization problem presents constraints on design variables. In particular, necessary conditions for optimal solutions are derived applying Pontryagin's Principle [18,19] and solved either analytically or numerically, depending on the complexity of the involved equations.
In this paper, for the sake of simplicity, we focus on the mechanical behavior of FG thick-walled cylinders. However, the considered optimization problems can be applied to FG rotating disks, shells of revolution, or other axisymmetric bodies. In the next section, the mechanical behavior of FG axisymmetric cylinders are briefly recalled. In particular, the governing equations referring to equilibrium, kinematic, and constitutive relations of linear plane elasticity theory are presented. Section 3 is dedicated to the formulation of the optimization problem, taking into account, in particular, the control input constraints and some possible cost functionals. In Section 4, a method to find a solution is described, containing some further discussion on the results. A special class of optimal solutions is thoroughly discussed in Section 5, highlighting optimal performances with respect to those obtained by means of classical property variations found in literature within a numerical example. Finally, conclusions are drawn in Section 6.

Governing Equations
Consider a uniformly pressurized cylinder of inner radius R i and outer radius R o and let p i be the internal pressure ( Figure 1a). The equilibrium equation of the infinitesimal element ( Figure 1b) in plane stress theory yields where σ r and σ θ denote radial and hoop stresses, both functions of r, while the derivative with respect to the radius is denoted with a "prime" symbol. The kinematic equations (or strain-displacement relations), which put into relation radial and circumferential strains ε r and ε θ with the radial displacement u, are Finally, the plane stress state linear constitutive relations (or Hooke's laws) which link the strains with stresses are where E is the Young's modulus, which is a function of r, ε y is the axial strain and ν is the Poisson coefficient, assumed constant along the radius for simplicity.
The radial strain in (2) can be written as which, together with (3), yields From (1), the hoop stress σ θ and its first derivative with respect to r are and respectively. Substituting (6) and (7) in (5) and rearranging the terms, one obtains the following (dimensionless) second-order differential equation for the radial stresŝ subject toσ (8) is commonly denoted as the Beltrami-Michell differential equation. Once solved, the hoop stress is computed using (6) and consequently the strains and the displacement from (3) and (2), respectively.

Formulation of the Optimal Control Problem
In the context of systems theory, a control problem consists in determining the values of the (input) control variables to move the state of the system from a given initial condition to a final one. More precisely, denoting by x and u the state and the control variables, respectively, and referring to the dynamics dx the control problem characterized by initial conditions x 0 , final conditions x T , admissible state space X and admissible control space V, consists in finding a control vector function v : The optimal control problem associated with a cost functional J(x, v) consists in finding, among all the solutions of the control problem, the one which minimizes J. Hence, in order to formulate an optimization problem, for the thick-walled cylinder, in the framework of dynamic programming, a state-space representation and a cost (or goal) functional are needed [18,19]. Note that, since in the governing differential Equation (8) derivatives are taken with respect to the (normalized) radial coordinater, the independent variable of the state-space representation is not time, as in Equation (10), butr. Therefore, instead of initial and final conditions, a set of boundary conditions at the internal and external surfaces of the wall are needed. Introducing the state variables x 1 =σ r , x 2 = dσ r /dr and x 3 = E/E i , the second order differential Equation (8) may be written as the first-order system where w is the control function, i.e., the derivative of Young's modulus with respect tor. This choice is justified since this latter appears explicitly in (8). Note that not all the boundary states are specified.
In particular, x 1 (1) and (9), yielding while x 2 (1) and x 2 (R o /R i ) are unknown and can be used to define the cost function, as shown below. Similarly, if the Young's modulus at the inner radius is fixed, then

Determination of the Cost Function
In dynamic programming theory and calculus of variations, cost functionals are mainly classified as Lagrange, Mayer, and Bolza functionals [20]: • A cost functional is of the Lagrange type if it consists in a distributed term associated with an integral over whole the considered interval. For system (11), it is associated with an integral over the interval [1, R i /R o ]. For instance, the minimization of the objective functional leads to minimizing the variance of the equivalent Tresca stresŝ keeping it as close as possible to the constant valueσ > 0 and making the stress as uniform as possible along the radius. • A cost functional is of the Mayer type if it consists in a function depending only on the final state conditions. Correspondingly, for system (11), a cost functional of such a kind depends only on the value of the state variables x 1 , x 2 , and x 3 at R o /R i . For instance, minimizing the functional leads to the minimization of the maximum equivalent stress for fixed a value of R i , under the hypothesis that the latter is strictly decreasing along the radius. • Finally, a cost functional of the Bolza type is the sum of a Mayer type and a Lagrange one. For instance, minimizing the functional leads to the minimization of the normalized displacement at the outer radius u(R o ) with respect to the internal pressure value, with respect to the inner radius and with respect to the Young's modulus at the inner boundary.
Given the boundary conditions (12) and (13), the solution of system (11) is unique for any choice of the input w; hence, the cost functionals J i (i = 1, 2, 3), although defined with respect to the state variables, are, in turn, functions of w.

Input Constraints
According to [21], there are a few optimization studies in which the manufacturability cost is taken into consideration. When considering not merely the theoretical possibility of designing FG materials but the technological processes employed to obtain them, it is reasonable to make the hypothesis that variations of a property, such as the Young's modulus along the thickness of a cylinder, are supposed to vary between two limit values. As a consequence, in the formulation of the optimal control problem, it is reasonable to assume the derivative of Young's modulus (with respect to the radius) w be constrained in an admissible range of values. More precisely, we assume, for all values of r, w ∈ [w − , w + ], where w − and w + can be deduced from fixed radial property variations or from technological process data.
The optimal control problem can now be stated formally. In the formulation of the problem, as well as in the computation of the solution, reference is made to the cost functional (14). However, solutions associated with the costs (15) or (16) can be obtained in a similar way. Problem 1. Given the dynamical system (11) and the boundary conditions (12) and (13), find the control function w : (14) is minimized.

Computation of the Solution
The following theoretical analysis is based on Berkovitz's framework [22]. Pontryagin's maximum principle applied to Problem 1 states that the optimal control function w, i.e., the one which minimizes the cost functional J 1 (w) is, among all admissible functions, the one which minimizes, for any value of r, the Hamiltonian function H defined by where x = (x 1 x 2 x 3 ) and p = (p 1 p 2 p 3 ) is the vector of the so-called co-state variables, all functions ofr. Note that the Hamiltonian function exhibits a linear dependence on the control function w, i.e., Equation (17) can be written as H(r, x, p, w) = a(r) + b(r)w where a(r) = (rx 2 −σ) 2 + p 1 (r)x 2 (r) − 3p 2 (r)x 2 (r) r and b(r) = p 3 (r) + p 2 (r) Since the problem is characterized by a closed control set, Pontryagin's Principle yields boundary solution for the minimization of (18). More precisely, the optimal control function w * is defined by which in the context of dynamic programming is usually called a bang-bang control, since the control function only assumes the extremal values of the admissible interval.
Since the control function is defined as the derivative of Young's modulus, optimal variations of E turn out to be, quite unexpectedly, piecewise linear with respect tor. The result is particularly interesting since the linearity of the Young's modulus is reasonably the easiest form of variation amenable for physical realization from the technological viewpoint.
Equation (21) does not yet provide the explicit expression of the optimal solution; in fact, it is clear that, in order to know the explicit value of w for any value ofr, one should know the value of b(r). In turn, the computation of b(r) through (20) requires the knowledge of the solution of the dynamical system (11) and of the differential equations for the co-state, which can be written explicitly as Boundary conditions for co-states are determined by the transversality conditions; by writing the cost functionals in the compact form a common expression of the transversality conditions is which, for the cost functional J 1 , read In the framework of control theory, the instants satisfying b = 0 are called switching instants since the control function switches instantaneously from an extreme value to the other. Borrowing this terminology, in the following, the values ofr satisfying b = 0 are called switching points.
The application of Pontryagin's Principle, therefore, leads to a system of six nonlinear and coupled first-order differential equations described by (11) and (23) that have to be solved by intervals for a constant value of w (in each interval either w ≡ w + or w ≡ w − ) and taking into account the six boundary conditions (12), (13) and (26).
To solve this boundary value problem (BVP), a numerical integration approach is needed. One way is to implement the so-called shooting method [24], i.e., to start from the inner radius and guess x 2 (1), p 1 (1), and p 3 (1), and then forward integrate (11) and (23) as an initial value problem (IVP) until the outer radius, where a check is made whether (12) and (26) are satisfied. If so, a solution is found, if not the initial guesses are adjusted. Other techniques are based on the so-called pseudospectral methods, where the optimal control problem is transcribed into a nonlinear constrained programming problem, whose solution can be computed by means of gradient methods. The numerical solution is out of the scope of the present study. However, the reader is addressed to [25,26], where a general background on numerical methods for optimal control problems is given.

The Case of a Single Switching Point
As mentioned above, information on switching points can be deduced from the number of roots of b. In fact, there are as many switches as the number of distinct (real) roots. In the following, we describe in detail the solution associated with the case of a single switching point. Moreover, we suppose that, in the formulation of Problem 1, the values of Young's moduli are fixed at both radii, i.e., that not only x 3 (1) is fixed but also x 3 (R o /R i ). This assumption allows one to define the linear variation, between the two boundary values, whose importance is twofold: (i) it will be used in a geometrical interpretation of admissible solutions and (ii) it will serve a comparison for the effectiveness of the bang-bang optimal solution. More precisely, letw denote the linear variation of Two situations may occur. Ifw / ∈ [w − , w + ], there is no admissible solution. In fact, as shown in Figure 2a, the line connecting E i with E o does not lie in either of the two angular sectors determined by w − (red dashed line) and w + (blue dashed line) and having the vertex in E i and E o , respectively. On the other hand, ifw ∈ [w − , w + ], the linear variation between E i and E o belongs to both the angular sectors (see Figure 2b) and is, therefore, an admissible solution; however, it is not optimal. Indeed, the optimal solution is one of the two solutions associated with the extremal values of w. Each of them is characterized by a different switching point and a different value of w in the two resulting subintervals; more precisely, they can be described as follows: • One solution (black solid line in Figure 2b) is characterized by a sub-interval (R i ,r 1 ) in which w * ≡ w + and a sub-interval (r 1 , R o ) in which w * ≡ w − . The switching pointr 1 is determined imposing In the following, we refer to this solution as the "P-M" solution.

•
The other solution (violet solid line in Figure 2b), which will be referred to as the "M-P" solution, is characterized by a sub-interval (R i ,r 2 ) in which w * ≡ w − and a sub-interval (r 2 , R o ) in which w * ≡ w + . The switching pointr 2 is determined imposing

A Numerical Example
To have an immediate insight of the importance of the result, we have simulated the application of the two bang-bang solutions to a particular instance of the equivalent stress minimization problem (functional J 2 ), and we have compared their performance with the ones exhibited by the simplest one, i.e., the linear variation, and one of the most used in the literature, i.e., the power law variation. We have considered a cylinder of internal radius R i = 20 mm and external radius from 30 mm to 50 mm (so that the relative external radius varies from 1.5 to 2.5). Internal and external pressures have been set to p i = 10 MPa and p o = 0 MPa, respectively. Different choices of FG materials are possible, depending on the manufacturing process to be used; nevertheless, in this case, internal and external values of the Young's modulus have been chosen to be E i = 2.1 × 10 5 MPa (steel) and E o = 3.9 × 10 5 MPa (alumina). This choice was justified by the fact that steel/alumina FG materials were adopted in several works dealing with FG material optimization, particularly in the case of axisymmetric structures (see: [6,27,28]). According to these values, the upper limit for w − is (see Equation (27)) = 4.5 × 10 3 MPa/mm while the lower limit for w + is If w − and w + do not fulfill these bounds, then we fall in the case represented in Figure 2a. In particular, we have considered three cases where w − = 4 × 10 3 MPa/mm, while w + = 20 × 10 3 MPa/mm, 24 × 10 3 MPa/mm and 28 × 10 3 MPa/mm, respectively. Radial and hoop stresses have been forecast by means of an axisymmetric FE analysis, whose details are omitted for brevity. However, the reader is addressed to [17], where similar computational aspects are reported. The variation of the maximum normalized Tresca stress σ T max /p i , as a function of the relative outer radius, is reported, for all the three cases and for both the bang-bang admissible variations, in Figure 3.
The P-M solution (dashed lines), which performs better than the M-P one (solid lines), is the optimal bang-bang solution. To evaluate its efficiency in reducing the stress, it is interesting to compare the stress values with the ones obtained for the linear (dash-dotted lines) and for the power law (dotted lines) variations of the Young's modulus. As one may note from Figure 3, the optimal bang-bang solution P-M still performs better than both the latter solutions.  To better characterize the optimal bang-bang solution, a numerical computation of the variation of the position of the switching points has been performed as a function of the relative thickness R o /R i , whose implicit expressions are (28) and (29). Results are reported, for the three cases, in  We can conclude that, at the cost of a little increase in complexity, the bang-bang solution-which is piecewise linear and, hence, rather simple and reasonably realizable-performs much better than the linear variation.

Conclusions
Material property variation in functionally graded materials has been reported in the context of dynamic programming framework. In particular, cylinders with variable Young's modulus along the radius subject to internal pressure have been considered. Several problem formulations have been stated and solutions are addressed by means of Pontryagin's Principle. Optimal Young's modulus turns out to be piecewise linear along the radius, the consequence of a bang-bang control scenario and is amenable for an easy physical realization.
The optimal solution requires one to resort to numerical methods. Comments on the switching points and possible extension including additional constraints are addressed. A numerical example has been performed to minimize the maximum equivalent stress within the assumption that a single switching point exists. The optimal solution performs much better than classic solutions present in literature.
Notwithstanding the idea is in its early formative stages, the authors strongly believe it brings novelty in the way of modeling property variation in functionally graded bodies as well as it can represent a touchstone for future developments and collaborations between mechanical designers and control theorists.