**1. Introduction**

Temperature effects on the plastic deformation have significant industrial and academic interests [1,2]. However, many studies in the literature make assumptions that the material properties are isotropic and temperature-independent (TI). For example, aluminum composite discs under thermal loading have been studied without using temperature-dependent (TD) material properties [3,4]. Considerations of elastic and plastic anisotropy are important when the deformation of textured metals or single crystals under thermal loading are in question [5]. Orthotropic plasticity material models have been extensively developed by Hill [6–8]. Anisotropic plasticity theory have been applied in many studies to understand directional dependent yielding phenomena. For example, Yoon10 et al. conducted research on the calibration of parameters used in anisotropic yield criterion from experimental tests in strongly textured aluminum sheets [9]. Numerical studies on predicting earing phenomena in anisotropic aluminum have been performed [10]. In addition, orthotropic plastic deformation in fiber-reinforced composite disc under spinning has been analyzed [11].

The importance of using temperature-dependent material properties in thermal loading analysis has been emphasized by Noda [12]. Thermomechanical responses of solid and hollow cylinders with the consideration of temperature-dependent material properties have been reported [13]. Realistic temperature functions to describe material properties at elevated temperatures for structural steels can be found in [14]. Elastoplastic stress analysis of thin discs with temperature-dependent material properties have been studied with analytical methods [15,16]. In addition, effects of thickness variations on the elastoplastic behavior of annular discs have been studied [17]. Although these studies consider the temperature-dependent material properties, they only deal with a system containing single material. Composite systems introduce additional complexity into the plasticity problem. Zarandi et al. examined the plastic responses of a composite disc, in two and three dimensions, under monotonic temperature loading with consideration of temperature dependent material properties [18]. In addition, the plasticity problem of a particular type of composite materials, termed functionally graded materials, under bending have recently been investigated [19]. The derived analytical solutions provide an efficient way in designing such materials.

In this work, we assume that the matrix material has orthogonal symmetry both in its elastic and plastic properties, such as a single crystal, metal with texture, or or fiber-reinforced composite materials. The matrix material is fully constrained on its outer rim, and its inner rim is in perfect bonding with an isotropic, purely elastic inclusion. The finite element method is adopted to conduct parametric studies on elastoplastic behavior and residual stress of the composite cylinder, analyzed under the plane-strain assumption, subjected to uniform, quasi-static thermal loading and unloading. Numerical schemes for solving elastoplasticity problems have been well established [20,21]. In this study, we conduct parametric studies to analyze residual stresses with software package [22]. Both temperature dependent and temperature independent material properties are considered. Effects of hardening and plane stress/strain are analyzed. Our numerical results may serve as reference data for experimental verifications or future analytical solutions to such a problem.

### **2. Theoretical and Numerical Considerations**

As shown in Figure 1a, the composite cylinder consists of an isotropic, purely elastic inclusion and elastoplastic matrix with the orthotropic symmetry both in its elastic and plastic behavior. The inclusion-matrix interface is assumed to be perfectly bonded. The outer boundary of the composite cylinder is fully clamped, and a uniform temperature field is quasi-statically applied to the composite cylinder. Figure 1b shows representative thermal loading cycles. The thermal loading parameter serves as loading steps in our analysis. In the figure, Δ *T* at Points B, D and F is in ratio of 1:1.25:1.5. The temperature differences at the three points are 700, 875 and 1050 ◦C, respectively. During the thermal loading or unloading cycles, residual stress fields may be developed at Points C, E and G. The physical properties of the isotropic elastic inclusion (Young's modulus *Ei* = 411 GPa, Poisson's ratio *νi* = 0.28, linear thermal expansion coefficient *αi* = 5.0 × 10−<sup>6</sup> <sup>K</sup>−1) are assumed to be temperature independent, but those of the elastoplastic matrix are temperature dependent. Physically it is envisioned that the inclusion is made of a ceramic material (0 < *r* < *a*), surrounded by metallic material (*a* < *r* < *b*), where *r* is the radial component of the polar coordinate system to describe the points in the domain. The elastoplastic problems are numerically analyzed via the finite element method in two dimensions.

**Figure 1.** (**a**) Schematic of the composite cylinder, with *a* = 0.3 m and *b* = 1 m, confined on the outer rim, *r* = 1 m, and (**b**) representative quasi-static temperature loading profile.

In the elastic region, the orthotropic constitutive relationship for the matrix is as follows.

$$
\varepsilon\_{ij} = S\_{ijkl}\sigma\_{kl} \qquad \text{or} \qquad \varepsilon\_{\mathfrak{m}} = S\_{\mathfrak{m}m}\sigma\_{\mathfrak{n}} \tag{1}
$$

where, at the reference temperature, *E*11 = 190, *E*22 = 200, *E*33 = 210, *ν*12 = 0.25, *ν*23 = 0.3, *ν*13 = 0.35, *G*12 = 75, *G*12 = 80, *G*12 = 85. Moduli are in units of GPa. Thermal expansion coefficient is assumed to have negligible orientational dependence, hence *αm* = 11.7 × 10−<sup>6</sup> in units of 1/K.

In the orthotropic symmetry, the Cartesian coordinates, (*<sup>x</sup>*1, *x*2, *<sup>x</sup>*3) or (*<sup>x</sup>*, *y*, *<sup>z</sup>*), are the rolling direction, the transverse direction and the normal direction, respectively. The Hill's orthogonal yield function is defined as

$$\mathcal{F}(\sigma\_{\overline{i}\overline{j}}) = \mathcal{F}\left(\sigma\_{22} - \sigma\_{33}\right)^2 + \mathcal{G}\left(\sigma\_{33} - \sigma\_{11}\right)^2 + H\left(\sigma\_{11} - \sigma\_{22}\right)^2 + 2L\sigma\_{23}^2 + 2M\sigma\_{13}^2 + 2N\sigma\_{12}^2 \tag{2}$$

It is assumed that at the reference temperature yield stresses are *σ*11 = 310, *σ*22 = 410, *σ*33 = 510, *σ*23 = 200, *σ*13 = 300, *σ*12 = 400 in units of MPa. The associate flow rule is as follows.

$$
\varepsilon\_{ij}^p = \lambda \frac{\partial \mathcal{F}}{\partial \sigma\_{ij}} \tag{3}
$$

Since the trance of the plastic strain tensor is zero due to the incompressibility assumption during plastic flow,

$$
\varepsilon\_{11}^p = \dots \left[ \text{ $ } \text{$  } G \left( \sigma\_{33} - \sigma\_{11} \right) + H \left( \sigma\_{11} - \sigma\_{22} \right) \right] \tag{4}
$$

$$
\varepsilon\_{22}^p = \varepsilon\_{21} \left[ F \left( \sigma\_{22} - \sigma\_{33} \right) - H \left( \sigma\_{11} - \sigma\_{22} \right) \right] \tag{5}
$$

$$
\varepsilon\_{33}^p = \,^2\Delta \left[ -F \left( \sigma\_{22} - \sigma\_{33} \right) + G \left( \sigma\_{33} - \sigma\_{11} \right) \right] \tag{6}
$$

 The consequence if the incompressibility assumption leads to,

$$
\varepsilon\_{11}^p + \varepsilon\_{22}^p + \varepsilon\_{33}^p = 0 \tag{7}
$$

The inter-relationships among Hill's parameters and directional yield stresses are as follows.

$$
\sigma\_{y23} = \sqrt{\frac{1}{2L}}, \qquad \sigma\_{y31} = \sqrt{\frac{1}{2M}}, \qquad \sigma\_{y12} = \sqrt{\frac{1}{2N}} \tag{8}
$$

and

$$F\_{\perp} = \frac{1}{2} \left( \frac{1}{\sigma\_{y22}^2} + \frac{1}{\sigma\_{y33}^2} - \frac{1}{\sigma\_{y11}^2} \right) \tag{9}$$

$$G\_{\parallel} = \frac{1}{2} \left( \frac{1}{\sigma\_{y33}^2} + \frac{1}{\sigma\_{y11}^2} - \frac{1}{\sigma\_{y22}^2} \right) \tag{10}$$

$$H\_{\parallel} = \frac{1}{2} \left( \frac{1}{\sigma\_{y11}^2} + \frac{1}{\sigma\_{y22}^2} - \frac{1}{\sigma\_{y33}^2} \right) \tag{11}$$

In this work, we specify the six yield stresses, instead the Hill's parameters. One may further define equivalent initial yield stress

$$
\sigma\_{y0} = \sqrt{\frac{3}{2(F+G+H)}}\tag{12}
$$

hence the Hill's effective stress

$$\sigma\_{\rm Hill}^2 = \sigma\_{y0}^2 \left[ \mathcal{F} \left( \sigma\_{22} - \sigma\_{33} \right)^2 + \mathcal{G} \left( \sigma\_{33} - \sigma\_{11} \right)^2 + H \left( \sigma\_{11} - \sigma\_{22} \right)^2 + 2L\sigma\_{23}^2 + 2M\sigma\_{13}^2 + 2N\sigma\_{12}^2 \right] \tag{13}$$

The plastic potential used in the isotropic hardening

$$\mathcal{F}\_h = \sigma\_{\text{Hill}} - \sigma\_y \tag{14}$$

where

$$
\sigma\_y = \sigma\_{y0} + \sigma\_h(\varepsilon\_{ep}) \tag{15}
$$

The hardening function *σh* depends on effective plastic strains *<sup>ε</sup>ep*. In the present analysis,

$$
\sigma\_h(\varepsilon\_{\varepsilon p}) = E\_{iso} \varepsilon\_{\varepsilon p} \tag{16}
$$

where

$$\frac{1}{E\_{\rm iso}} = \frac{1}{E\_{\rm This}} - \frac{1}{E} \tag{17}$$

and *ETiso* is the isotropic tangent modulus and *E* the effective Young's modulus if the material is elastically anisotropic. In this work, we set *ETiso* = 20 MPa if linear hardening is considered. The local effective plastic strain is defined as follows.

$$
\dot{\epsilon}\_{ep} = \sqrt{\frac{2}{3} \dot{\epsilon}\_{ij}^p \dot{\epsilon}\_{ij}^p} \tag{18}
$$

The von Mises effective stress is defined as follows in terms of deviatoric stress tensor *sij*, or its second invariant *J*2 [5].

$$
\sigma\_{\text{misces}} = \sqrt{3j\_2(s\_{ij})} = \sqrt{\frac{3}{2}s\_{ij}s\_{ij}} \quad \text{and} \quad s\_{ij} = \sigma\_{ij} - \frac{1}{3}\sigma\_{kk}\delta\_{ij} \tag{19}
$$

here *δij* is the Kronecker delta function, and the Einstein summation rule for the indices is applied. The temperature functions, shown in Equations (20)–(23), are used in this study. They are similar to those in Argeso and Eraslan ([13]), but slightly modified.

$$f\_{\sigma}(T) = \sigma\_y(T) / \sigma\_0 = 1 + T / [600 \times \ln(T / 1630)] \tag{20}$$

$$f\_E(T) = E(T) / E\_0 = 1 + T / \left[ 2000 \times \ln(T / 1800) \right] \tag{21}$$

$$f\_{\nu}(T) = \nu(T)/\nu\_{0} = 1 + 2.5 \times 10^{-4} T - 2.5 \times 10^{-7} T^{2} \tag{22}$$

$$f\_a(T) = a(T)/a\_0 = 1 + 2.56 \times 10^{-4} T - 2.14 \times 10^{-7} T^2 \tag{23}$$

where the applied temperature difference *T* is in units of ◦C and reference temperature is 22 ◦C. A graphical representation of the temperature functions is shown in Figure 2. For the case that the material is elastically and plastically isotropic, at the reference temperature, the Young's modulus, yield stress, Poisson's ratio and linear thermal expansion coefficient of the matrix material are assumed to be *E*0 = 200 GPa, *σ*0 = 410 MPa, *ν*0 = 0.3, and *α*0 = 11.7 × 10−<sup>6</sup> per ◦C, respectively. When the material is orthotropic, the above-mentioned temperature functions are applied to the material parameters at the reference temperature, listed after Equation (1). In other words, temperature-dependent yield stress *σ*(*T*) = *σ*0 *fσ*(*T*), Young's modulus *E*(*T*) = *E*0 *fE*(*T*), Poisson's ratio *ν*(*T*) = *ν*0 *fν*(*T*) and linear thermal expansion coefficient *α*(*T*) = *α*0 *fα*(*T*). This choice of material parameters is representative for studying the temperature dependent elastoplastic materials in general. The deformation process is assumed to be quasi-static throughout this work.

**Figure 2.** (**a**) Temperature functions for the material properties and (**b**) representative finite element mesh used in this study.

A representative finite element mesh is shown in Figure 2b with the inclusion radius *a* = 0.3 m and the radius to the outer boundary of the cylinder is *b* = 1 m. The number of two-dimensional quadratic serendipity elements used in the analysis was about 6000, and the number of degrees of freedom (d.o.f.) was about 170,000, including the internal d.o.f. for plasticity. We adopted COMSOL ([22]) software for the finite element calculations.

### **3. Results and Discussion**

### *3.1. Residual Stress Analysis Under Elastic-Perfectly Plastic Assumption*

When the plane-plane composite cylinder, with elastic-perfectly plastic model, under temperature cyclic loading, Figure 3 shows the residual stress at *r* = 0.35 m, near the inclusion-matrix boundary, in the TD case under the maximum applied temperature differences Δ *T* = 270, 337.5, 405 ◦C. All material parameters are assumed to be temperature dependent. The maximum applied temperature differences are labeled as the letters 'B', 'D' and 'F' in the Figure 1, as schematics. Since the chosen loading magnitudes are large, reversed plasticity occurs during unloading. The stress magnitudes at the polar angle *θ* = 45◦ direction are larger than those at *θ* = 0◦ due to the chosen plasticity parameters in the Hill's model having larger yield stress when *θ* = 45◦.

Under the TI assumption, Figure 4 shows the residual stress in the orthotropic cylinder under maximum temperature differences Δ*T* = 270, 337.5, 405 ◦C, which are the same as those used in the TD case. Since the loading magnitudes are kept the same in Figures 3 and 4, direct comparisons to exhibit the the differences between TD and TI assumptions can be accomplished. In general, since TI does not reduce yield stress at high temperature, it predicts higher stresses. Furthermore, it can be seen that the residual stresses, at the loading parameters about 2.5, 4.5 and 7, are mildly developed during reversed plasticity in the TI case since at high temperatures yield stress is not deduced. These loading parameters respectively correspond to the letters 'C', 'E' and 'G', shown in the schematics in Figure 1b. However, in the TD case, yield stress is largely reduced at high temperatures as shown in Figure 2a. This large reduction in yield stress causes strong reversed plasticity in the unloading phases. Hence, the TD residual stresses in Figure 3, at the loading parameters about 2.5, 4.5 and 7, are much larger than those in the TI case.

**Figure 3.** Von Mises residual stress at *r* = 0.35 m in the temperature-dependent (TD) case with the polar angle (**a**) *θ* = 0◦ and (**b**) *θ* = 45◦. Dashed lines indicate the applied temperature difference.

**Figure 4.** Von Mises residual stress at *r* = 0.35 m in the temperature-independent (TI) case with the polar angle (**a**) *θ* = 0◦ and (**b**) *θ* = 45◦. Dashed lines indicate the applied temperature difference.

### *3.2. Effects of Selective Temperature Dependent Material Properties on Residual Stress*

In the previous section, the results are obtained with all material properties being temperature dependent. In this section, we examine the effects of selective TD material properties on residual stress. Linear hardening with *ET* = 20 MPa is assumed and its temperature dependence is assumed to be

*fE*(*T*). Figure 5a shows the residual stress in the orthotropic cylinder after the maximum temperature difference Δ*T* = 400 ◦C loading with all material parameters being temperature dependent, i.e., the four temperature functions listed in Equations (20)–(23) are used in the analysis. When only yield stress is considered to be temperature dependent, the corresponding residual stress is shown in Figure 5b. As can be seen their von Mises residual stress distributions are similar, but their magnitudes are distinct. Furthermore, the range for the residual stress 'plateau' is shorter when all material properties are temperature dependent. The 'plateau' is slightly inclined due to the small linear hardening *ET*. As for comparisons, the residual stress distribution for the TI case is shown in Figure 5c. Due to no reduction in yield stress in the TI case as temperature increases, sharp residual stress distribution is developed near the inclusion-matrix interface. It is remarked that residual stresses are self-equilibrated inside the material. However, the von Mises residual stresses do not show this trend since they have averaged according to Equation (19).

**Figure 5.** Residual stress distribution for (**a**) all material parameters being temperature dependent, (**b**) only yield stress being temperature dependent and (**c**) all material parameters being temperature independent.

### *3.3. Residual Stress Analysis with Linear Hardening*

When the plane-strain composite cylinder under maximum temperature differences Δ*T* = 400, 500, 600 ◦C, the residual stresses, at *r* = 0.35 m, in each loading cycle are shown in Figure 6. Both plastic anisotropy and isotropy are compared. Linear hardening is assumed, as in previous section. For the TD case, all material properties are assumed to be temperature dependent in this section. It can be seen that the TI case predicts large residual stress due to reversed plasticity in the unloading phase. Due to orthotropic elastic constants, the numerical values between the isotropic and orthotropic case can only be compared qualitatively.

**Figure 6.** Residual stress near the inclusion-matrix interface at *r* = 0.35 for the (**a**) TD and (**b**) TI case.

To observe the orientational dependence of residual stress in the cylinder, Figure 7 shows the stress contour in the cylinder after unloading from maximum temperature difference Δ*T* =400 ◦C. Color bars indicate the von Mises residual stress in units of MPa. The TI case shows weaker developments in residual stress due to yield stress not being reduced at elevated temperature. Due to the matrix being both elastically and plastically orthotropic, the isotropic inclusion also shows Under the same loading/unloading conditions, Figure 8 shows the residual stress developed in the isotropic composite cylinder. As expected, no orientational dependence is observed in the residual stress distribution.

**Figure 7.** Von Mises residual stress in the orthotropic composite cylinder with the (**a**) TD and (**b**) TI material properties.

**Figure 8.** Von Mises residual stress in the isotropic composite cylinder with the (**a**) TD and (**b**) TI material properties.

In order to better observe the directional dependence of the residual stress, shown in Figure 7, we plot the residual stress at *r* = 0.35 m, near the inclusion-matrix interface indicated by a black circle in the figure, around the circumferential direction, as shown in Figure 9. The circumferential direction is indicated by the polar angle whose zero value is at the x-axis. It is remarked that for isotropic plasticity such a plot would show horizontal lines only; no orientational dependence. As can be seen in Figure 9, there are significant differences in the von Mises residual stress between the TD and TI case. In the TD case, due to significant reduction in yield stress, the reversed plasticity during unloading phase induces complex residual stress distribution.

**Figure 9.** Von Mises residual stress at *r* = 0.35 m around the polar angle in the orthotropic cylinder for the (**a**) TD and (**b**) TI case.

Although many practical problems need to be analyzed in three space dimensions, long cylinders can be analyzed in two dimensions under the plane strain assumption. Through the basic research on the parametric analysis, it is demonstrated that the purpose of this study is to examine the residual stress in the orthotropic cylinder under thermal loading with the consideration of temperature dependent material properties. If materials have strong temperature dependence, one needs to be aware significant differences between the TD and TI case. In the future, considerations of temperature rates are required in order to realistically model the material responses.

Possible real-world applications of the present analysis are machine parts with anisotropic characteristics under repeated temperature loading/unloading cycles or textured alloys during metal forming processes. In addition, fiber-reinforced composite materials with fibers being arranged as concentric rings may require the consideration of both elastic and plastic orthotropy, as studied in the present work. This type of composite materials may be used in machinery, civil engineering or biomedical engineering.

The development of residual stresses in materials depend on a variety of variables [23]. If the residual stresses are to be minimized for certain applications, suitable selection of materials and boundary conditions, as well as the methodology of analyzing the residual stresses under given conditions, need to be considered. From the present study, to avoid complex residual stress developments in the TD case under thermal loading, one may consider to choose a material with less anisotropy and temperature sensitivity. For the TI case, one may either choose more isotropic material or use multiple layers near the interface to reduce the orientational dependence in residual stress. The multiple layers with suitable design and material selections may reduce the stress in the matrix, hence reduce yielding and consequent the developments of residual stresses after unloading. The multiple layer method may also work for the TD case. Further analysis is required.
