**1. Introduction**

Thin rotating disks are used in many applications such as energy storage devices; gyroscopic control devices for ships, submarines, aircrafts, rockets, and missiles; high-speed gears; and turbine rotors [1]. Moreover, rotational autofrettage has been recently proposed [2] as a new technique for producing compressive residual stresses. A purely elastic solution for isotropic rotating disks has been given [3], and a comprehensive overview of the problem of an elastic rotating disk up through to the late 1990s has been provided [4]. Further, an elastic solution for arbitrarily functionally graded polar orthotropic rotating disks has been recently proposed [5].

A grea<sup>t</sup> number of elastic/plastic solutions have been also reported in the literature. In most cases it has been assumed that plastic yielding is controlled by the Tresca yield criterion. A review of such solutions has been provided [6]. A few solutions for the deformation theory of plasticity based on the von Mises yield criterion are also available; a review of these solutions has been given [7]. In the case of the flow theory of plasticity, the finite difference method has usually been adopted for determining the distribution of strains [8–10]. An efficient method that advances the analytical treatment of elastic/plastic rotating disks has been proposed [11] for the von Mises yield criterion and its associated flow rule. However, the general idea of the method can be extended to other yield criteria with no difficulty. In particular, it has been demonstrated [12] where a model of pressure-dependent plasticity has been adopted.

Solutions for anisotropic materials are of special importance because even mild plastic anisotropy has an amplified effect on residual stress and strain distributions under certain conditions [13]. The proposed orthotropic yield criterion [14] is often used to model plastic anisotropy in rotating disks, for example [15,16]. An important type of anisotropy is polar anisotropy (see, for example, [16–20]). Therefore, the present paper deals with polar orthotropic disks obeying the yield criterion [14]. The method developed in [11] is employed. This allows a semi-analytic solution to be found. Comparison with the finite difference solution presented in [10] is made. It is shown that the finite difference solution is not accurate enough for some applications. This demonstrates an advantage of using the method [11] as compared with the finite difference method.

### **2. Statement of the Problem**

A detailed description of the boundary value problem under consideration can be found in many works (see, for example, [11]). The boundary value problem is solved in a cylindrical coordinate system (*r*, *θ*, *z*) whose *z* axis coincides with the axis of symmetry of a thin annular rotating disk of constant thickness. The assumption of constant thickness is often adopted in theoretical analyses of rotating disks [21–28]. The outer and inner radii of the disk are denoted as *b*0 and *a*0, respectively. The angular velocity of the disk is *ω*. The boundary value problem is illustrated in Figure 1. Strains are small. The disk has no stress at *ω* = 0. The normal stresses in the cylindrical coordinate system, *σr*, *σθ*, and *σz*, are the principal stresses. The solution of the boundary value problem is independent of the polar angle. The state of stress is plane, *σz* = 0. The component of the acceleration vector in the circumferential direction is neglected. The stress boundary conditions are:

$$
\sigma\_{\mathbb{T}} = 0 \tag{1}
$$

for *r* = *a*0 and *r* = *b*0. In general, the disk consists of two regions: elastic and plastic. Hooke's law connects the elastic strains and stresses. In particular,

$$
\varepsilon\_r^\varepsilon = \frac{\sigma\_\mathcal{r} - \nu \sigma\_\theta}{E}, \quad \varepsilon\_\theta^\varepsilon = \frac{\sigma\_\theta - \nu \sigma\_\mathcal{r}}{E}, \quad \varepsilon\_z^\varepsilon = -\frac{\nu (\sigma\_\mathcal{r} + \sigma\_\theta)}{E}. \tag{2}
$$

Here, *ν* is Poisson's ratio and *E* is Young's modulus, *εr* is the radial strain, *εθ* is the circumferential strain, and *εz* is the axial strain. The superscript e denotes the elastic part of the strain and will denote the elastic part of the strain rate. The orthotropic yield criterion proposed in reference [14] and its associated flow rule are adopted in the plastic region. It is assumed that the principal axes of anisotropy coincide with coordinate curves of the cylindrical coordinate system. Under plane stress conditions, the yield criterion adopted reads (*G* + *<sup>H</sup>*)*σ*2*r* + (*F* + *<sup>H</sup>*)*σ*2*θ* − 2*Hσrσθ* = 1 where *G*, *H*, and *F* are anisotropic constants. It is convenient to rewrite this criterion as:

$$
\sigma\_r^2 + \frac{\sigma\_\theta^2}{\eta\_1^2} - \frac{\eta \sigma\_r \sigma\_\theta}{\eta\_1} = \sigma\_0^2 \tag{3}
$$

where

$$\eta = \frac{2H}{\sqrt{(G+H)(H+F)}}, \quad \eta\_1 = \frac{\sqrt{G+H}}{\sqrt{H+F}}, \quad \eta\_2 = \frac{2F}{(H+F)\sqrt{4-\eta^2}}, \quad \sigma\_0 = \frac{1}{G+H}.\tag{4}$$

**Figure 1.** Schematic diagram of a rotating annular disk.

The yield criterion (Equation (3)) reduces to the von Mises yield criterion at *F* = *G* = *H*. In this case, *σ*0 is the tensile yield stress. Let .*εpr* , .*εpθ* , and .*εpz* be the plastic strain rates. In the case under consideration, the associated flow rule can be written as [29]:

$$\dot{\varepsilon}\_r^p = \lambda \left( \sigma\_r - \frac{\eta}{2\eta\_1} \sigma\_\theta \right), \quad \dot{\varepsilon}\_\theta^p = \frac{\lambda}{\eta\_1} \left( \frac{\sigma\_\theta}{\eta\_1} - \frac{\eta}{2} \sigma\_r \right), \quad \dot{\varepsilon}\_z^p = -\lambda \left[ \left( 1 - \frac{\eta}{2\eta\_1} \right) \sigma\_r + \left( \frac{1}{\eta\_1} - \frac{\eta}{2} \right) \frac{\sigma\_\theta}{\eta\_1} \right] \tag{5}$$

where *λ* is a non-negative multiplier. The superimposed dot denotes the time derivative at fixed *r* and the superscript *p* denotes the plastic part of the strain rate and will denote the plastic part of the strain. The total strains and strain rates in the plastic region are:

$$\begin{array}{ll}\varepsilon\_r = \varepsilon\_r^\varepsilon + \varepsilon\_{r\prime}^p & \varepsilon\_\theta = \varepsilon\_\theta^\varepsilon + \varepsilon\_{\theta\prime}^p & \varepsilon\_z = \varepsilon\_z^\varepsilon + \varepsilon\_{z\prime}^p\\ \dot{\varepsilon}\_r = \dot{\varepsilon}\_r^\varepsilon + \dot{\varepsilon}\_r^p & \dot{\varepsilon}\_\theta = \dot{\varepsilon}\_\theta^\varepsilon + \dot{\varepsilon}\_\theta^\theta & \dot{\varepsilon}\_z = \dot{\varepsilon}\_z^\varepsilon + \dot{\varepsilon}\_z^p \end{array} \tag{6}$$

The equilibrium equation is of the form:

$$\frac{\partial \sigma\_r}{\partial r} + \frac{\sigma\_r - \sigma\_\theta}{r} = -\varsigma \omega^2 r. \tag{7}$$

Here, *ς* is the mass density of the material.

The boundary value problem is classified by the following dimensionless quantities:

$$
\rho = \frac{r}{b\_0}, \quad \Omega = \frac{\varsigma \omega^2 b\_0^2}{\sigma\_0}, \quad a = \frac{a\_0}{b\_0}, \quad k = \frac{\sigma\_0}{E}. \tag{8}
$$

Since the material model is rate-independent, the time derivative can be replaced with the derivative with respect to Ω. The derivatives of strain components with respect to Ω are denoted as:

$$\begin{split} \xi^{r}\_{r} &= \frac{\partial \varepsilon\_{r}}{\partial \Omega}, \quad \xi^{\theta}\_{\theta} = \frac{\partial \varepsilon\_{\theta}}{\partial \Omega}, \quad \xi^{z}\_{z} = \frac{\partial \varepsilon\_{z}}{\partial \Omega}, \\\xi^{x}\_{r} &= \frac{\partial \varepsilon^{r}\_{r}}{\partial \Omega}, \quad \xi^{x}\_{\theta} = \frac{\partial \varepsilon^{\theta}\_{\theta}}{\partial \Omega}, \quad \xi^{x}\_{z} = \frac{\partial \varepsilon^{x}\_{z}}{\partial \Omega}, \\\xi^{p}\_{r} &= \frac{\partial \varepsilon^{p}\_{r}}{\partial \Omega}, \quad \xi^{p}\_{\theta} = \frac{\partial \varepsilon^{p}\_{\theta}}{\partial \Omega}, \quad \xi^{p}\_{z} = \frac{\partial \varepsilon^{p}\_{z}}{\partial \Omega}. \end{split} \tag{9}$$

Then, the equation of strain rate compatibility is equivalent to:

$$
\rho \frac{\partial \zeta\_{\theta}}{\partial \rho} = \zeta\_r - \zeta\_{\theta}. \tag{10}
$$

Using Equation (8), Equation (7) can be rewritten as:

$$\frac{\partial \sigma\_r}{\sigma\_0 \partial \rho} + \frac{\sigma\_r - \sigma\_\theta}{\sigma\_0 \rho} = -\Omega \rho. \tag{11}$$

### **3. Solution at Loading**

### *3.1. Purely Elastic Solution*

The general purely elastic solution for a rotating disk is given, for example, in [3]. Using Equation (8), the solution satisfying the boundary condition (Equation (1)) at *ρ* = 1 is represented as:

$$\begin{split} \frac{\sigma\_{\tau}}{\sigma\_{0}} &= B \left( \frac{1}{\rho^{2}} - 1 \right) + \frac{\Omega(3+\nu)}{8} (1 - \rho^{2}), \\ \frac{\sigma\_{0}}{\sigma\_{0}} &= -B \left( \frac{1}{\rho^{2}} + 1 \right) + \frac{\Omega(1+3\nu)}{8} \left( \frac{3+\nu}{1+3\nu} - \rho^{2} \right), \\ \frac{\tau\_{\tau}}{k} &= \frac{8B \left[ 1 + \nu - (1-\nu)\rho^{2} \right] + \Omega(1-\nu) \left[ 3 + \nu - 3(1+\nu)\rho^{2} \right] \rho^{2}}{8\rho^{2}}, \\ \frac{\tau\_{\theta}}{k} &= \frac{-8B \left[ 1 + \nu + (1-\nu)\rho^{2} \right] + \Omega(1-\nu) \left[ 3 + \nu - (1+\nu)\rho^{2} \right] \rho^{2}}{8\rho^{2}}, \\ \frac{\tau\_{\tau}}{k} &= \frac{\nu}{4} \left\{ 8B - \Omega \left[ 3 + \nu - 2(1+\nu)\rho^{2} \right] \right\}. \end{split} \tag{12}$$

Here, *B* is a constant of integration. Using the boundary condition (Equation (1)) at *ρ* = *a*, it is possible to find that:

$$B = -\frac{\Omega(3+\nu)a^2}{8}.\tag{13}$$

Eliminating *B* in Equation (12) by means of Equation (13) provides the solution for the stresses and strains in the purely elastic disk. Let Ω*e* be the value of Ω at which the plastic region starts to develop from the surface *ρ* = *a*. Substituting the values of *σr* and *σθ* at *ρ* = *a* into the yield condition Equation (3) and using Equation (4) yields:

$$
\Omega\_{\varepsilon} = \frac{4\eta\_1}{3 + \nu + a^2(1 - \nu)}.\tag{14}
$$

In what follows, it is assumed that Ω > Ω*<sup>e</sup>*.

The solution (Equation (12)) is valid in the elastic region of the elastic/plastic disk. However, Equation (13) is not valid in this case because one of the boundary conditions in Equation (1) should be replaced with the condition that the radial and circumferential stresses are continuous across the elastic/plastic boundary.

### *3.2. Elastic/Plastic Stress Solution*

The elastic/plastic stress solution is available [30]. For completeness, this solution is outlined below. The yield condition (Equation (3)) is satisfied by the following substitution:

$$
\sigma\_{\mathbb{T}}/\sigma\_0 = 2\cos\psi/\sqrt{4-\eta^2}, \quad \sigma\_{\mathbb{R}}/\sigma\_0 = \left(\eta\cos\psi/\sqrt{4-\eta^2} + \sin\psi\right)\eta\_1 \tag{15}
$$

where *ψ* is a new function of *ρ* and Ω. Substituting Equation (15) into Equation (11) and using Equation (4) leads to:

$$\frac{2\sin\psi}{\sqrt{4-\eta^2}}\frac{\partial\psi}{\partial\rho} + \frac{(\eta\_2\cos\psi - \eta\_1\sin\psi)}{\rho} = \Omega\rho. \tag{16}$$

The boundary condition to this equation is determined from the boundary condition (Equation (1)) at *ρ* = *a* and Equation (15) as:

$$
\psi = \frac{\pi}{2} \tag{17}
$$

for *ρ* = *a*. It has been taken into account here that *σθ* > 0 at *ρ* = *a*. Let *ρc* be the dimensionless radius of the elastic/plastic boundary. The value of *ψ* at *ρ* = *ρc* is denoted by *ψ<sup>c</sup>*. The continuity of the radial and circumferential stresses across the elastic/plastic boundary, together with Equations (12) and (15), leads to: 

$$\begin{split} B\left(\frac{1}{\rho\_c^2} - 1\right) + \frac{\Omega(3+\nu)}{8} (1 - \rho\_c^2) &= 2\cos\psi\_c / \sqrt{4 - \eta^2}, \\ -B\left(\frac{1}{\rho\_c^2} + 1\right) + \frac{\Omega(1+3\nu)}{8} (\frac{3+\nu}{1+3\nu} - \rho\_c^2) &= \left(\eta\cos\psi\_c / \sqrt{4 - \eta^2} + \sin\psi\_c\right)\eta\_1. \end{split} \tag{18}$$

Eliminating *B* between these equations yields:

$$\begin{split} & \frac{2\cos\psi\_c}{\sqrt{4-\eta^2}} + \frac{\left(\eta\cos\psi\_c / \sqrt{4-\eta^2} + \sin\psi\_c\right)\left(1-\rho\_c^2\right)\eta\mathbf{1}}{\left(1+\rho\_c^2\right)} \\ & -\frac{\Omega\left(1-\rho\_c^2\right)}{8} \left[\frac{\left(1+3\nu\right)}{\left(1+\rho\_c^2\right)}\left(\frac{3+\nu}{1+3\nu} - \rho\_c^2\right) + 3 + \nu\right] = 0. \end{split} \tag{19}$$

For a given value of Ω, Equation (19) and the solution of Equation (16) supply the system of equations for *ρc* and *ψ<sup>c</sup>*. Then, *B* can be determined from Equation (18). The stress distribution in the elastic region, *ρc* ≤ *ρ* ≤ 1, follows from Equation (12). The stress distribution in the plastic region, *a* ≤ *ρ* ≤ *ρ<sup>c</sup>*, is readily determined from Equation (15) and the solution of Equation (16). The latter is in parametric form with *ψ* being the parameter. The plastic region occupies the entire disk when *ρc* = 1. Let Ω *p* be the corresponding value of Ω. This value can be found numerically using the dependence of *ρc* on Ω known from Equation (19) and the solution of Equation (16). In particular, it is seen from Equation (19) that *ψc* = *π*/2 at *ρc* = 1. This condition and Equation (17) allow the value of Ω *p* and a constant of integration to be determined from the solution of Equation (16).

### *3.3. Elastic/Plastic Strain Solution*

Eliminating *B* in Equation (12) by means of the solution of Equations (18) and (19) supplies the strain solution in the elastic region as follows in terms of Ω and *ρ*. Replacing the plastic strain rates in Equation (5) with the corresponding quantities introduced in Equation (9) and eliminating *λ* between the resulting equations leads to:

$$\xi\_r^p = \xi\_\theta^p \frac{(2\eta\_1\sigma\_r - \eta\sigma\_\theta)\eta\_1}{(2\sigma\_\theta - \eta\eta\_1\sigma\_r)}, \quad \xi\_z^p = \xi\_\theta^p \frac{[(2\eta\_1 - \eta)\eta\_1\sigma\_r + (2 - \eta\eta\_1)\sigma\_\theta]}{(\eta\eta\_1\sigma\_r - 2\sigma\_\theta)}.$$

The stresses in these equations can be expressed in terms of *ψ* by means of Equation (15). Then,

$$\xi\_{\tau}^{\mathbb{Z}^p} = \xi\_{\theta}^p \frac{\left(\sqrt{4-\eta^2}\cos\psi - \eta\sin\psi\right)\eta\_1}{2\sin\psi}, \quad \xi\_{\tau}^p = -\xi\_{\theta}^p \frac{\left[\eta\_1\sqrt{4-\eta^2}\cos\psi + (2-\eta\eta\_1)\sin\psi\right]}{2\sin\psi}.\tag{20}$$

The elastic portion of the strain components in the plastic region is determined from Equations (2), (8) and (15) as:

$$\begin{aligned} \frac{c\_r^\varepsilon}{k} &= \frac{(2 - \eta \eta \chi \nu) \cos \psi}{\sqrt{4 - \eta^2}} - \eta\_1 \nu \sin \psi, \quad \frac{c\_\theta^\varepsilon}{k} = \frac{(\eta \eta\_1 - 2 \nu) \cos \psi}{\sqrt{4 - \eta^2}} + \eta\_1 \sin \psi, \\\frac{c\_\theta^\varepsilon}{k} &= \frac{\nu (\eta \eta\_1 + 2) \cos \psi}{\sqrt{4 - \eta^2}} + \nu \eta\_1 \sin \psi. \end{aligned} \tag{21}$$

Differentiating these expressions with respect to Ω and using Equation (9) yields:

$$\begin{split} \frac{\partial \xi}{\partial \xi} &= -\left[ \frac{(2-\eta\eta\_1\nu)\sin\psi}{\sqrt{4-\eta^2}} + \eta\_1\nu\cos\psi \right] \frac{\partial \psi}{\partial \Gamma \Omega}, \quad \frac{\xi\_\theta^x}{k} = \left[ \eta\_1\cos\psi - \frac{(\eta\eta\_1 - 2\nu)\sin\psi}{\sqrt{4-\eta^2}} \right] \frac{\partial \psi}{\partial \Gamma \Omega}, \\\frac{\xi\_\xi^x}{k} &= \nu \left[ \eta\_1\cos\psi - \frac{(\eta\eta\_1 + 2)\sin\psi}{\sqrt{4-\eta^2}} \right] \frac{\partial \psi}{\partial \Gamma \Omega}. \end{split} \tag{22}$$

Taking into account Equation (6), Equation (10) can be rewritten as *ρ∂ξθ*/*∂ρ* = *ξ p r* + *ξe r* − *ξ p θ* − *ξe θ*. Eliminating *ξ p r* in this equation by means of Equation (20) results in:

$$\rho \frac{\partial \xi\_{\theta}}{\partial \rho} = \frac{\left[\eta\_1 \sqrt{4 - \eta^2} \cos \psi - (\eta \eta\_1 + 2) \sin \psi\right]}{2 \sin \psi} \xi\_{\theta}^p + \xi\_r^\kappa - \xi\_{\theta}^\kappa.$$

This equation and Equation (6) combine to give the following equation for *ξθ*:

$$\rho \frac{\partial \underline{\boldsymbol{\varepsilon}}\_{\theta}}{\partial \rho} = \frac{\left[\eta\_1 \sqrt{4 - \eta^2} \cos \psi - (\eta \eta\_1 + 2) \sin \psi\right]}{2 \sin \psi} \underline{\boldsymbol{\varepsilon}}\_{\theta} + \frac{\eta\_1 \left[\eta \sin \psi - \sqrt{4 - \eta^2} \cos \psi\right]}{2 \sin \psi} \underline{\boldsymbol{\varepsilon}}\_{\nu}^{\varepsilon} - \underline{\boldsymbol{\varepsilon}}\_{\theta}^{\varepsilon} \tag{23}$$

where *ξe r* and *ξe θ* should be eliminated by means of Equation (22). It is therefore evident that Equation (23) is a linear ordinary differential equation. In order to solve this equation, it is necessary to express the derivative *∂ψ*/*∂*Ω involved in Equation (22) in terms of *ψ* or/and *ρ*. Following the method proposed in [8], Equation (16) is differentiated with respect to Ω. Then,

$$\frac{2\sin\psi}{\sqrt{4-\eta^2}}\frac{\partial\chi}{\partial\rho} + \left(\frac{2\cos\psi}{\sqrt{4-\eta^2}}\frac{\partial\psi}{\partial\rho} + \frac{\eta\_2\sin\psi + \eta\_1\cos\psi}{\rho}\right)\chi - \rho = 0\tag{24}$$

where *χ* = *∂ψ*/*∂*Ω. Using Equation (16), the derivative *∂ψ*/*∂ρ* involved in Equation (24) can be expressed in terms of *ψ* and *ρ*. Then, Equation (24) becomes:

$$\frac{2\sin\psi}{\sqrt{4-\eta^2}}\frac{\partial\chi}{\partial\rho} + \left\{\frac{\eta\_2\sin\psi + \eta\_1\cos\psi}{\rho} + 2\cos\psi \left[\Omega\rho - \frac{(\eta\_2\cos\psi - \eta\_1\sin\psi)}{2\rho\sin\psi}\right] \right\}\chi - \rho = 0. \tag{25}$$

It follows from the boundary condition (Equation (17)) that *∂ψ*/*∂*Ω = 0 at *ρ* = *a*. Therefore,

$$
\chi = 0\tag{26}
$$

at *ρ* = *a*. This is the boundary condition to Equation (25), which is a linear ordinary differential equation for *χ*. This equation should be solved numerically. Once *χ* has been found from Equation (25), it is possible to determine *ξe r* and *ξe θ* involved in Equation (23) as functions of *ψ* and *ρ* by means of Equation (22). Having *ψ* as a function of *ρ* due to the solution of Equation (16), it is possible to represent the coefficients of Equation (23) as functions of *ρ*. Then, this ordinary differential equation can be solved numerically with no difficulty. The boundary condition to Equation (23) follows from the continuity of the circumferential strain rate and, therefore, *ξθ* across the elastic/plastic boundary *ρ* = *ρ<sup>c</sup>*. The value of *ξθ* on the elastic side of this boundary is determined from Equation (12) as:

$$\frac{\rho\_c^2}{k} = \frac{(1-\nu)\left[3+\nu-(1+\nu)\rho\_c^2\right]}{8} - \frac{\left[1+\nu+(1-\nu)\rho\_c^2\right]}{\rho\_c^2}\frac{dB}{d\Omega}.\tag{27}$$

Then, the boundary condition to Equation (23) is

$$
\xi\_{\theta} = \xi\_{c} \tag{28}
$$

for *ρ* = *ρ<sup>c</sup>*. Once the solution of Equation (23) has been found, the total circumferential strain in the plastic region is determined by integration of *ξθ* with respect to Ω at a given point *ρ* = *ρ<sup>t</sup>*. The maximum value of Ω is denoted as <sup>Ω</sup>*f* and the corresponding value of *ρc* as *ρf* . It is assumed that <sup>Ω</sup>*f* < <sup>Ω</sup>*p*. It is obvious that *a* ≤ *ρt* < *ρf* . The value of Ω and the value of the circumferential strain on the elastic side of the elastic/plastic boundary at *ρc* = *ρt* are denoted as Ω*t* and *Et θ*, respectively. Then, it follows from the definition for *ξθ* given in Equation (9) that:

$$\varepsilon\_{\theta} = \int\_{\Omega\_{l}}^{\Omega\_{f}} \xi\_{\theta} d\Omega + E\_{\theta}^{t} \,. \tag{29}$$

The value of *Etθ*is determined from Equation (12) as:

$$\frac{E\_{\theta}^{t}}{k} = \frac{-8B\_{t}\left[1+\nu+(1-\nu)\rho\_{t}^{2}\right]+\Omega\_{t}(1-\nu)\left[3+\nu-(1+\nu)\rho\_{t}^{2}\right]\rho\_{t}^{2}}{8\rho\_{t}^{2}}.\tag{30}$$

Here *Bt* is the value of *B* at *ρc* = *ρ<sup>t</sup>*. This value follows from Equations (18) and (19). Equations (6) and (22) combine to give:

$$\xi\_{\theta}^{p} = \xi\_{\theta} - k \left[ \eta\_{1} \cos \psi - \frac{(\eta \eta\_{1} - 2\nu) \sin \psi}{\sqrt{4 - \eta^{2}}} \right] \chi \,. \tag{31}$$

Substituting Equation (31) into Equation (20) leads to:

$$\begin{split} \frac{\xi\_{\mathbb{Z}}^{\mathbb{P}}}{k} &= \frac{\eta\_{1}}{2} \left\{ \frac{\tilde{\varepsilon}\_{\theta}}{k} - \left[ \eta\_{1} \cos \psi - \frac{(\eta \eta\_{1} - 2\upsilon) \sin \psi}{\sqrt{4 - \eta^{2}}} \right] \chi \right\} \frac{\left( \sqrt{4 - \eta^{2}} \cos \psi - \eta \sin \psi \right)}{\sin \psi}, \\\ \frac{\xi\_{\mathbb{Z}}^{\mathbb{Z}}}{k} &= -\frac{1}{2} \left\{ \frac{\tilde{\varepsilon}\_{\theta}}{k} - \left[ \eta\_{1} \cos \psi - \frac{(\eta \eta\_{1} - 2\upsilon) \sin \psi}{\sqrt{4 - \eta^{2}}} \right] \chi \right\} \frac{\left[ \eta\_{1} \sqrt{4 - \eta^{2}} \cos \psi + (2 - \eta \eta\_{1}) \sin \psi \right]}{\sin \psi}. \end{split} \tag{32}$$

Hence,

$$\frac{\varepsilon\_r^p}{k} = \int\_{\Omega\_t}^{\Omega\_f} \frac{\chi\_r^p}{k} d\Omega\_r \quad \frac{\varepsilon\_z^p}{k} = \int\_{\Omega\_t}^{\Omega\_f} \frac{\chi\_z^p}{k} d\Omega \;. \tag{33}$$

Here, the integrands are known functions of Ω. In particular, *ξ pr* and *ξ pz* are first eliminated by means of Equation (32). Then, the solutions of Equations (16), (23) and (25) are used to represent *ψ*, *ξθ*, and *χ*, respectively, as functions of Ω at any value of *ρ* = *ρ<sup>t</sup>*. Therefore, the integrals involved in Equation (33) can be evaluated numerically. The total radial and axial strains in the plastic zone are found by means of Equation (6) where the plastic parts are given by Equation (33) and the elastic parts by Equation (21).

### **4. Distribution of Residual Stresses and Strains**

It is assumed that unloading is purely elastic (i.e., the distribution of residual stresses found by means of Hooke's law for the increments of stresses does not violate the yield criterion in the range *a* ≤ *ρ* ≤ 1). This assumption should be verified a posteriori. The residual stresses are determined as:

$$
\sigma\_r^{\rm res} = \sigma\_r^f + \Delta \sigma\_r, \quad \sigma\_\theta^{\rm res} = \sigma\_\theta^f + \Delta \sigma\_\theta \,. \tag{34}
$$

Here, *σfr* and *σfθ* are the radial and circumferential stresses, respectively, at the end of loading. These stresses were found in the previous section. Δ*σr* and Δ*σθ* are the increments of the radial and circumferential stresses, respectively, in the course of the process of unloading. Analogously, the residual strains are determined as:

$$
\varepsilon\_r^{\rm res} = \varepsilon\_r^f + \Delta \varepsilon\_r \quad \varepsilon\_\theta^{\rm res} = \varepsilon\_\theta^f + \Delta \varepsilon\_\theta \quad \varepsilon\_z^{\rm res} = \varepsilon\_z^f + \Delta \varepsilon\_z. \tag{35}
$$

Here, *εfr*, *εfθ* , and *εfz* are the total radial, circumferential, and axial strains, respectively, at the end of loading. These strains were found in the previous section. Δ*εr*, Δ*εθ*, and Δ*εz* are the increments of the radial, circumferential, and axial strains, respectively, in the course of the process of unloading. Since the process of unloading is assumed to be purely elastic, the increments of the strains are related by Hooke's law to the increments of the stresses:

$$
\Delta \varepsilon\_r^{\varepsilon} = \frac{\Delta \sigma\_r - \nu \Delta \sigma\_\theta}{E}, \quad \Delta \varepsilon\_\theta^{\varepsilon} = \frac{\Delta \sigma\_\theta - \nu \Delta \sigma\_r}{E}, \quad \Delta \varepsilon\_z^{\varepsilon} = -\frac{\nu \left(\Delta \sigma\_r + \Delta \sigma\_\theta\right)}{E}. \tag{36}
$$

Therefore, the solution (Equations (12) and (13)) in which Ω should be replaced with −Ω*<sup>f</sup>* is valid. As a result,

$$\begin{aligned} \frac{\Delta r\_{\!
\!
}}{\tau\_{0}} &= B\_{f} \left( \frac{1}{\rho^{2}} - 1 \right) - \frac{\Omega\_{f} (3 + \nu)}{8} (1 - \rho^{2}),\\ \frac{\Delta \sigma\_{\!
\!
}}{\tau\_{0}} &= -B\_{f} \left( \frac{1}{\rho^{2}} + 1 \right) - \frac{\Omega\_{f} (1 + 3\nu)}{8} \left( \frac{3 + \nu}{1 + 3\nu} - \rho^{2} \right),\\ \frac{\Delta \sigma\_{\!
\!
}}{\tilde{k}} &= \frac{8B\_{f} \left[ 1 + \nu - (1 - \nu)\rho^{2} \right] - \Omega\_{f} (1 - \nu) \left[ 3 + \nu - 3(1 + \nu)\rho^{2} \right] \rho^{2}}{8\rho^{2}},\\ \frac{\Delta \epsilon\_{\!
\!
}}{\tilde{k}} &= \frac{-8B\_{f} \left[ 1 + \nu + (1 - \nu)\rho^{2} \right] - \Omega\_{f} (1 - \nu) \left[ 3 + \nu - (1 + \nu)\rho^{2} \right] \rho^{2}}{8\rho^{2}},\\ \frac{\Delta \epsilon\_{\!
\!
}}{\tilde{k}} &= \frac{\nu}{4} \left\{ 8B\_{f} + \Omega\_{f} \left[ 3 + \nu - 2(1 + \nu)\rho^{2} \right] \right\} \end{aligned} \tag{37}$$

where

$$B\_f = \frac{\Omega\_f (3+\nu) a^2}{8}.\tag{38}$$

Substituting the solution found in the previous section together with Equations (37) and (38) into Equation (35) yields the distribution of residual stresses and strains. It follows from Equation (3) that the yield criterion is not violated after unloading if:

$$
\left(\frac{\sigma\_r^{\rm res}}{\sigma\_0}\right)^2 + \left(\frac{\sigma\_\theta^{\rm res}}{\sigma\_0 \eta\_1}\right)^2 - \frac{\eta}{\eta\_1} \left(\frac{\sigma\_r^{\rm res}}{\sigma\_0}\right) \left(\frac{\sigma\_\theta^{\rm res}}{\sigma\_0}\right) - 1 \le 0. \tag{39}
$$

Since the distribution of the residual stresses has been found, this inequality can be verified with no difficulty.

### **5. Illustrative Example**

Equations (16), (23) and (24) were solved numerically for several materials whose anisotropic coefficients, given in [31,32], are shown in Table 1. It is assumed that *ν* = 0.3 in all cases. The value of k introduced in Equation (8) is immaterial. In particular, assume that the solution for a disk of given geometry and physical properties has been found. Then, simple scaling of this solution supplies the solutions for similar disks of material with the same Poisson's ratio and anisotropic coefficients but any value of k. The numerical solution is illustrated for an *a* = 0.5 disk at <sup>Ω</sup>*f* = 1.7. The distribution of stresses is depicted in Figures 2 and 3 and residual stresses in Figures 4 and 5. The associated total and plastic strain fields are shown in Figures 6–11. The variation of the residual strains with *ρ* is depicted in Figures 12–14. The solution for the residual stresses was used in conjunction with Equation (39) to verify that the process of unloading is purely elastic. It is evident from Figure 3 that the effect of plastic anisotropy on the radius of the elastic plastic boundary is very significant and, as a result, so is the effect of plastic anisotropy on the distribution of stresses, strains, residual stresses, and residual strains. On the other hand, the yield loci for the anisotropic parameters considered are depicted in Figure 15. It is seen from this figure that the yield loci for the materials considered are not very different (except for the yield locus for AA3104). This sensitivity of solutions to the yield locus requires very accurate numerical methods for calculating stress and strain fields. Therefore, it is of interest to compare the present solution and a finite difference solution. The present solution involves fewer approximations than finite difference solutions because the derivative *∂ψ*/*∂*Ω is found from Equation (24) without any discretization with respect to Ω (i.e., with respect to time). Therefore, it is natural to assume that the present solution is more accurate. In [10], several finite difference solutions were found for an *a* = 3/7 disk assuming that *ρf* = 5/7 (in our nomenclature). Using the anisotropic constants adopted

in [10], the distribution of the circumferential strain was determined using the method proposed in the present paper. A comparison of the magnitude of circumferential strain at *ρ* = *a* and *ρf* = 5/7 (in our nomenclature) predicted by the new method and that found in [10] is presented in Table 2. It can be seen from this table that the accuracy of the finite difference solution may be insufficient for practical applications. In particular, Δ shown in Table 2 is defined as:

$$
\Delta = \frac{|\varepsilon\_{\theta, FDM} - \varepsilon\_{\theta, N}|}{\varepsilon\_{\theta, N}} \times 100\%
$$

where *εθ*,*FDM* is the total circumferential strain at *ρ* = *a* and *ρf* = 5/7 found in [10] and *εθ*,*<sup>N</sup>* is the total circumferential strain at *ρ* = *a* and *ρf* = 5/7 found in the present paper.


**Table 1.** Anisotropic coefficients of several materials.

**Table 2.** Comparison of the total circumferential strain at *ρ* = *a* and *ρf* = 5/7 found by the two methods.


**Figure 2.** Effect of anisotropic properties on the distribution of the radial stress in an *a* = 0.5 disk.

**Figure 3.** Effect of anisotropic properties on the distribution of the circumferential stress in an *a* = 0.5 disk.

**Figure 4.** Effect of anisotropic properties on the distribution of the residual radial stress in an *a* = 0.5 disk.

**Figure 5.** Effect of anisotropic properties on the distribution of the residual circumferential stress in an *a* = 0.5 disk.

**Figure 6.** Effect of anisotropic properties on the distribution of the radial strain in an *a* = 0.5 disk.

**Figure 7.** Effect of anisotropic properties on the distribution of the circumferential strain in an *a* = 0.5 disk.

**Figure 8.** Effect of anisotropic properties on the distribution of the axial strain in an *a* = 0.5 disk.

**Figure 9.** Effect of anisotropic properties on the distribution of the radial plastic strain in an *a* = 0.5 disk.

**Figure 10.** Effect of anisotropic properties on the distribution of the circumferential plastic strain in an *a* = 0.5 disk.

**Figure 11.** Effect of anisotropic properties on the distribution of the axial plastic strain in an *a* = 0.5 disk.

**Figure 12.** Effect of anisotropic properties on the distribution of the residual radial strain in an *a* = 0.5 disk.

**Figure 13.** Effect of anisotropic properties on the distribution of the residual circumferential strain in an *a* = 0.5 disk.

**Figure 14.** Effect of anisotropic properties on the distribution of the residual axial strain in an *a* = 0.5 disk.

**Figure 15.** Yield loci for several materials.

## **6. Conclusions**

A semi-analytic solution for the stresses and strains within a rotating elastic/plastic polar orthotropic annular disk was found under plane stress. The range of validity of plane stress solutions was determined in [33] by comparing such solutions with 3-D finite element solutions. The distribution of residual stresses and strains was determined as well. The constitutive equations consist of Hooke's law, the orthotropic yield criterion proposed in [14], and the associated flow rule. Therefore, the equations to be solved involve the strain rate tensor. This greatly adds to the difficulties of the solution as compared with the constitutive equations that relate the stresses to the strains (or allow for the strains to be immediately found by integrating relations between strain rate components). In order to facilitate numerical solution, the method developed in [11] was adopted. As a result, it is only necessary to use numerical methods for solving ordinary differential equations and to evaluate ordinary integrals, even though the solution depends on two independent variables—Ω and *ρ*.

It is seen from Equations (12), (21), (22), (32) and (33) that the parameter k introduced in Equation (8) is immaterial in the sense that scaling of any strain solution for a disk of given radius, Poisson's ratio, and anisotropic parameters provides the solutions for similar disks of material with the same Poisson's ratio and anisotropic parameters but any value of k.

It is known that numerical codes should be verified before their use in applications [34–36]. The present solution is useful for this purpose since ordinary differential equations can be solved numerically with a very high accuracy with no difficulty. In particular, comparison with the finite difference solution proposed in [10] was made. It was shown that the accuracy of the finite difference solution for the total circumferential strain at the inner radius of the disk may be insufficient (Table 2).

The solution found is for a rate-independent model of plasticity. However, in many cases, solutions for rate-dependent plasticity are required [37]. Moreover, disks of varying thickness and disks made of functionally graded materials are widely used in industry [38]. The general approach used in the present paper can be extended to at least some of these cases. This will be the subject of a subsequent investigation.

**Author Contributions:** All three authors participated in the research and in the writing of this paper.

**Funding:** This research received no external funding. **Acknowledgments:** This work was performed while W. Jeong was a visiting research professor at Beihang University, Beijing, China.

**Conflicts of Interest:** The authors declare no conflict of interest.
