Understanding the thermal problem of variable gradient functionally graded plate based on hybrid numerical method under linear heat source

Abstract

The thermal problem of functionally graded materials (FGM) under linear heat source is studied by a hybrid numerical method. The accuracy of the analytical method and the efficiency of the finite element method are taken into account. The volume fraction of FGM in the thickness direction can be changed by changing the gradient parameters. Based on the weighted residual method, the heat conduction equation under the third boundary condition is established. The temperature distribution of FGM under the action of linear heat source is obtained by Fourier transform. The results show that the closer to the heat source it is, the greater the influence of the heat source is and the influence of the heat source is local. The temperature change trend of the observation points is consistent with the heat source, showing a linear change. The results also show that the higher the value of gradient parameter is, the higher the temperature of location point is. The temperature distribution of observation points is positively correlated with gradient parameter. When the gradient parameter value exceeds a certain value, it has a little effect on the temperature change in the model and the heat conduction in the model tends to be pure metal heat conduction, the optimal gradient parameters combined the thermal insulation property of ceramics and the high strength toughness of metals are obtained.

Introduction

Functionally graded material (FGM) ¹ is a new kind of non-uniform composite material. In contrast to traditional composite materials, there is no sudden change in the material performance. FGM is widely used in the thermal domain since the structure can relieve thermal stress mutation and avoid thermal stress concentrations.

In recent years, researchers have studied the properties of FGM in analytical method. Sankar and Tzeng ² studied the thermal stress of FGM and derived an analytical solution of the model temperature distribution and thermal stress distribution. Alibeigloo ³ assumed that the physical parameters of functionally graded structures varied exponentially along the thickness direction, analytical solutions of the stress and temperature fields of the structure were derived. Shao ⁴ studied the axisymmetric static thermoelastic problem of a simply supported FGM cylinder using the method of differential equation sequence. Xu et al. ⁵ used the plane steady-state heat conduction equation and the method of separating variables, the analytical solutions of the steady-state temperature field of FGM plate at a range of temperatures were derived. Since the nonlinear of FGM, the analytical methods are often difficult to obtain. Therefore, the numerical methods^6,7 have been widely applied. Kim ⁸ presented graded finite elements within the framework of a generalized isoparametric formulation. A multiple isoparametric finite element method had been used by Zou et al. ⁹ In recent years, more numerical methods ¹⁰ have been proposed. Sladeka et al. ¹¹ used the local boundary integral equation method to study the transient heat conduction of FGM. Long ¹² put forward a basic concept and numerical implementation of a local Petrov-Galerkin method for solving the elasticity problem. Paulino et al. ¹³ and Sutradhar et al. ¹⁴ introduced the Galerkin boundary element method into the calculation of FGM steady-state temperature field, and then solved the FGM transient heat conduction problem by combining Laplace transform with the Galerkin boundary element method. Later, they solved the transient heat conduction problem of FGM using the boundary element method. ¹⁵ Mirzaei and Dehghan ¹⁶ used the moving least squares (MLS) to approximate in both time and space domains and then studied the heat transfer efficiency. Hassanzadeh-Aghdam and Ansari ¹⁷ studied the influence of agglomeration of SiC nanoparticles on the bulk thermal conducting behaviour of aluminium (Al)-based nanocomposites by using a micromechanics-based hierarchical technique. Ansari and Hassanzadeh-Aghdam ¹⁸ studied elastic, thermoelastic and viscoelastic properties of functionally graded carbon nanotube reinforced polymer nanocomposites using a 3-dimensional micromechanics.

The analytical solution is very difficult to be obtained, and its accuracy is reduced by using numerical methods. The hybrid numerical method (HNM) is a semi-analytical and semi-numerical method. For the thermal problems of functionally graded materials, the finite element method is used in one direction and Fourier transform is used in the other direction. It combines the accuracy of the analytical method with the efficiency of numerical method. This work has followed earlier studies that have shown the heat conduction of functionally graded structures with H(t) and exp(t).^19,20 heat source using hybrid numerical method, and its correctness has been verified by Liu and Xi ²¹ but the linear heat source is not mentioned. The purpose of this work is to explore the heat conduction of FGM with variable gradient parameters under the linear heat source. Based on the HNM, the influences of variable gradient parameters for a functionally graded plate on temperature distribution under linear heat source are obtained.

The establishment of the heat conduction theory of HNM under linear heat source

The length of the FGM plate in the $x$ and $y$ direction is assumed infinite, and the total thickness in the $z$ direction is $H$ . The thermal physical property in the $z$ direction is assumed to be continuous. The FGM plate is divided into $N$ elements along $z$ direction. The material in the plane along the thickness direction is homogeneous and isotropic. Each element has a thickness of $h$ and is divided into three nodal lines, namely the upper, middle and lower computational nodal line, and there is one computational temperature on each nodal line. When no heat source is applied, the temperature of FGM plate is zero. Then a linear heat source $T_f$ is applied at $x=0$ on the upper surface of the plate and the heat transfer coefficient between the upper surface and surrounding medium is $\alpha$ . The heat conduction problem is simplified to a two-dimensional problem that is related to only $x$ and $z$ , and the analysis model is shown in Figure 1.

Figure 1.

Heat conduction analysis model of the FGM plate.

The temperature field of FGM plate is given as ¹⁹

$$T=\mathbf{N}\mathbf{\Phi }$$

(1)

where T is the temperature function, $\mathbf{N}$ is a quadratic polynomial interpolation function in $z$ direction as following

The thickness direction is discretized by interpolation, and the interpolation shape function in the thickness direction is derived by the construction principle of displacement function.

$$\mathbf{N}=\left[\begin{array}{ccc}(1-3\overline{z}+2{\overline{z}}^2)\hfill & 4(\overline{z}-{\overline{z}}^2)\hfill & (2{\overline{z}}^2-\overline{z})\hfill \end{array}\right]$$

(2)

in which $\overline{z}=z/h$ , $h$ is the element thickness.

$\mathbf{\Phi }$ is a matrix of nodal temperature, which are functions of the coordinate $x$ and time $t$ at $z=0$ , $z=0.5h$ and $z=h$ , denoted as follows

$$\mathbf{\Phi }={({\varphi }_L,{\varphi }_M,{\varphi }_U)}^{\mathrm{T}}$$

(3)

Where ${\varphi }_L,{\varphi }_M,{\varphi }_U$ are respectively the $\mathbf{\Phi }$ value at the lower, middle and upper node surface.

The test function of the weighted residual value method is set as N^T, then the weighted residual value equation composed of boundary conditions and heat conduction equation is expressed as

$$\begin{array}{l}\underset{0}{\overset{h}{\int }}{\mathbf{N}}^{\mathrm{T}}\left[\eta \rho \frac{\partial T}{\partial t}-\left[{\gamma }_0\frac{{\partial }^{2}T}{\partial x^2}+{\gamma }_0\frac{{\partial }^{2}T}{\partial z^2}\right]-q_v\right]dz\\ +{{\gamma }_0\frac{\partial T}{\partial z}|}_{\mathrm{\Gamma }}+{\alpha (T-T_f)|}_{\mathrm{\Gamma }}=0\end{array}$$

(4)

here $\eta$ , $\rho$ , ${\gamma }_0$ , $q_v$ , $\alpha$ , $T$ and $t$ are material specific heat, mass density, thermal conductivity, internal heat source, heat transfer coefficient, computational temperature and calculation time, respectively. $T_f$ is the heat source temperature, subscript $\mathrm{\Gamma }$ denotes the value at the boundary.

The initial temperature is zero, it can be denoted as

$${T|}_{t=0}=0$$

(5)

equation(4)is simplified as

$$\mathbf{C}\stackrel{·}{\mathbf{\Phi }}+{\mathbf{K}}_D\mathbf{\Phi }=\mathbf{V}$$

(6)

Substitution of equation (1) into equation (4), obtain

$$\mathbf{C}={\int }_0^h{\mathbf{N}}^{\mathrm{T}}\eta \rho \mathbf{N}\mathrm{dz}=\frac{\eta \rho h}{30}\left[\begin{array}{ccc}4\hfill & 2\hfill & -1\hfill \\ \multicolumn{1}{c}{2}& 16\hfill & 2\hfill \\ \multicolumn{1}{c}{-1}& 2\hfill & 4\hfill \end{array}\right]$$

(7)

$${\mathbf{K}}_D=-{\mathbf{A}}_1·\frac{{\partial }^2}{\partial x^2}+{\mathbf{A}}_0$$

(8)

$${\mathbf{A}}_1=\frac{h{\gamma }_0}{30}\left[\begin{array}{ccc}4\hfill & 2\hfill & -1\hfill \\ \multicolumn{1}{c}{2}& 16\hfill & 2\hfill \\ \multicolumn{1}{c}{-1}& 2\hfill & 4\hfill \end{array}\right]$$

(9)

$${\mathbf{A}}_0=\alpha \left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -1\hfill \end{array}\right]+\frac{{\gamma }_0}{3h}\left[\begin{array}{ccc}7\hfill & -8\hfill & 1\hfill \\ \multicolumn{1}{c}{-8}& 16\hfill & -8\hfill \\ \multicolumn{1}{c}{1}& -8\hfill & 7\hfill \end{array}\right]$$

(10)

$$\mathbf{V}=\frac{q_{v}h}{6}\left[\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{4}\\ \multicolumn{1}{c}{1}\end{array}\right]+{\alpha T_f|}_{\mathrm{\Gamma }}$$

(11)

where the subscript $D$ denotes the matrix ${\mathbf{K}}_D$ is a differential operator matrix.

Assembled element equation to whole equation and applied the Fourier transform to equation (6), By using Fourier transform, the time domain problem is transformed into the wave number domain problem, and the differential term about x in the differential operator matrix is eliminated, then a set of system equations is as follows

$${\mathbf{C}}_g{\stackrel{·}{\stackrel{~}{\mathbf{\Phi }}}}_g+{\mathbf{K}}_{\mathrm{Dg}}{\tilde{\mathbf{\Phi }}}_g={\tilde{\mathbf{V}}}_g$$

(12)

Matrix ${\mathbf{K}}_{\mathrm{Dg}}$ is given by

$${\mathbf{K}}_{\mathrm{Dg}}={{\zeta }_x}^2{\mathbf{A}}_{1g}+{\mathbf{A}}_{0g}.$$

(13)

which is a constant matrix for given ${\zeta }_x$ . ${\stackrel{.}{\stackrel{~}{\mathbf{\Phi }}}}_g,{\tilde{\mathbf{\Phi }}}_g,{\tilde{\mathbf{V}}}_g$ are the Fourier transform of ${\stackrel{·}{\mathbf{\Phi }}}_g,{\mathbf{\Phi }}_g,{\mathbf{V}}_g$ , respectively. The subscript g denotes the matrix is whole matrix in computation field. Matrices ${\mathbf{C}}_g,{\mathbf{A}}_{1g},{\mathbf{V}}_g$ are obtained by assembling $\mathbf{C},{\mathbf{A}}_1,\mathbf{V}$ from all the elements. Since there are three nodal lines in one element, the dimension of ${\mathbf{C}}_g,{\mathbf{A}}_{1g}$ matrices is $M\times M(M=2n+1)$ , and the dimension of ${\mathbf{V}}_g$ matrices is $M\times 1(M=2n+1)$ , where $n$ is the number of the layer elements.

The temperature field in the space-time domain is obtained by inverse Fourier transform.

Solve equation (12), the result can be obtained as following.

$${\tilde{\mathbf{\Phi }}}_g=\sum _{m=1}^M{e}^{-{\varpi }_{m}t}(\int \frac{{\mathbf{ѱ}}_m^L{\tilde{\mathbf{V}}}_g}{{\mathbf{ѱ}}_m^L{\mathbf{C}}_g{\mathbf{ѱ}}_m^R}{e}^{{\varpi }_m\mathrm{dt}}\mathrm{dt}+O){\mathbf{ѱ}}_m^R$$

(14)

in which ${\varpi }_m$ and ${\mathbf{ѱ}}_m$ is the eigenvalue and the eigenvector of homogeneous equation of equation (12), respectively the superscripts $L$ and $R$ represent the left eigenvector and right eigenvector, respectively, $O$ is a constant.

For a given linear temperature heat source

$$T_f=T_{0f}t$$

(15)

where $T_{0f}$ is the coefficient factor of the heat source.

Regardless of the volume heat source, ${\mathbf{V}}_g$ can be written in terms of equation (11)

$${\mathbf{V}}_g={\alpha T_f|}_{\mathrm{\Gamma }g}$$

(16)

The Fourier transform is applied to equation (16), then

$${\tilde{\mathbf{V}}}_g={\alpha {\tilde{T}}_{0f}t|}_{\mathrm{\Gamma }g}$$

(17)

Substitution of equation (17) into equation (14), and using initial condition, the result of ${\tilde{\mathbf{\Phi }}}_g$ is solved as

$${\tilde{\mathbf{\Phi }}}_g=\sum _{m=1}^M\frac{{{{\mathbf{ѱ}}_m}^L\alpha {\tilde{T}}_{0f}|}_{\mathrm{\Gamma }g}({\varpi }_{m}t-1+e^{-{\varpi }_{m}t})}{{{\mathbf{ѱ}}_m}^L{\mathbf{C}}_g{{\mathbf{ѱ}}_m}^R{{\varpi }_m}^2}{{\mathbf{ѱ}}_m}^R$$

(18)

Performing the inverse Fourier transform to equation (18) and substituting into equation (1), the temperature field of FGM plate in linear heat source can be obtained.

Heat conduction under linear heat source

Parameter settings

The FGM is considered to be synthesized using a combination of a ceramic (thermal insulation) and metal (thermal conductor). The upper surface of the FGM is considered to be a pure ceramic material, namely zirconia ( $\mathrm{Zr}{\mathrm{O}}_2$ ), while the lower surface is a pure metal material ( $\mathrm{Ti}-6\mathrm{Al}-4\mathrm{V}$ ), the properties of intermediate material involve gradient changes in the thickness direction. The material properties of thermal conductivity, specific heat and mass density of $\mathrm{Zr}{\mathrm{O}}_2$ and $\mathrm{Ti}-6\mathrm{Al}-4\mathrm{V}$ are shown in Table 1.

Table 1.

Material properties of $\mathrm{Zr}{\mathrm{O}}_2$ and $\mathrm{Ti}-6\mathrm{Al}-4\mathrm{V}$ .

Material	Thermal conductivity $\gamma$ ( $\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1}$ )	Specific heat $\eta$ ( $\mathrm{J}·\mathrm{k}{\mathrm{g}}^{-1}·{\mathrm{K}}^{-1}$ )	Mass density $\rho$ ( $\mathrm{kg}·{\mathrm{m}}^{-3}$ )
$\mathrm{Zr}{\mathrm{O}}_2$	2.09	456.7	5331
$\mathrm{Ti}-6\mathrm{Al}-4\mathrm{V}$	7.50	537.0	4420

The temperature variation of the linear heat source over time can be obtained from Eq. (15), the value of $T_{0f}$ is taken as 330K, as shown in Figure 2. The time domain and space domain are separately discretized. The computing time is divided into 200 steps in 2 s, and the space domain is discretized into 1000 points for consideration of 0 to 10 H space domain in the x direction. The thickness of the $z$ direction is $H=0.06\mathrm{m}$ and divided into 40 elements. The heat transfer coefficient between the upper surface and surrounding medium is $\alpha =5*10^{-12}\mathrm{W}·{\mathrm{m}}^{-2}·{\mathrm{K}}^{-1}$

Figure 2.

Temperature variation of the linear heat source over time.

Establishment of variable gradient parameter model

Assume that the volume fraction of ceramic matrix material along z direction in FGM mathematical model is as follows.

$$V_C={(1-a(1-\frac{z}{H})+b{(1-\frac{z}{H})}^c)}^p$$

(19)

where $a,b,c,p$ are gradient parameters.

The gradient parameters for $a=1,b=0,p=(0,+\infty )$ are studied. The change of the volume fraction in the model is controlled by the gradient parameter p, and different volume fraction can be obtained by changing the p value. It can be seen that when the p value is about equal to 0, the V_C value is about equal to 1 and the material is close to the pure ceramic material. When the p value tends to infinity, the V_C value is about 0 and the material is close to the pure metal material. The volume fractions of the ceramics under different gradient parameters in the thickness are shown in Figure 3.

Figure 3.

The volume fraction of ceramic under different gradient parameters in thickness direction.

The main physical parameters considered in this paper are thermal conductivity, specific heat and mass density of FGM. These depend on the volume fraction of the ceramics in the model. The change trend of physical parameters is consistent. Here, the thermal conductivity change with different gradient parameters is shown as an example.

Assume that the thermal conductivity distribution in z direction is as follows.

$$\gamma (z)=({\gamma }_{\mathrm{C}}-{\gamma }_{\mathrm{M}})V_{\mathrm{C}}+{\gamma }_{\mathrm{M}}$$

(20)

where ${\gamma }_{\mathrm{C}}$ is the thermal conductivity of $\mathrm{Zr}{\mathrm{O}}_2$ , ${\gamma }_{\mathrm{M}}$ is the thermal conductivity of $\mathrm{Ti}-6\mathrm{Al}-4\mathrm{V}$ , and $V_C$ is the volume fraction of $\mathrm{Zr}{\mathrm{O}}_2$ at the relative position.

The thermal conductivity changes in thickness direction of FGM under different gradient parameters is shown in Figure 4.

Figure 4.

Thermal conductivity of FGM under different gradient parameters in thickness direction.

The result discussion of the heat conduction

Three observation position (a_i, b_i, c_i, i=1, 41) are investigated. They are 0.2H, 0.4H and 1.7H away from the heat source in the x direction, respectively. Different surface has different property, as can be shown in Figure 5.

Figure 5.

Schematic diagram of the observation position a_i, b_i, c_i of the FGM.

The temperature distribution of the upper surface and the lower surface at 1.0, 1.5 and 2.0 s is shown in Figures 6 and 7 while the gradient parameter is 5.0. It can be seen that the longer the heat source acts, the higher the temperature is. The temperature of 0 to 2 H distance from the heat source in x direction quickly drops, the temperature change of 2 to 10 H distance from the heat source in x direction is very small, and it fluctuates around lower values and then tends to zero. This indicates that the heat source has a great influence on the temperature distribution near the location and a very limited influence on the farther location. The influence of the heat source on FGM is localised, which is similar to the influence of the force field described by Saint Venant theory on a structure. It also can be seen that the temperature change trend of the upper surface and the lower surface is the same, the temperature of the lower surface is far lower than that of the upper surface and the upper surface plays a good thermal insulation effect.

Figure 6.

The temperature distribution of the upper surface at different times (p = 5).

Figure 7.

The temperature distribution of the lower surface at different times (p = 5).

The temperature change of observation points a₃₁ and a₁₁ under different gradient parameters is shown in Figures 8 and 9. It can be seen that the change trend of the temperature at different observation points is consistent with the heat source temperature function, showing a linear change. Figure 8 shows that when the gradient parameter values are 0.5, 2.0 and 5.0, the peak temperature of the observation point is 73.12, 192.83 and 211.91 K, the larger the value of the gradient parameter is, the higher the peak temperature of the observation point is and the temperature distribution has a positive correlation with the value of the gradient parameter. This is because the smaller the gradient parameter is, the larger the volume fraction of the ceramic is, which greatly reduces the thermal conductivity of the composite. Taking a₃₁ as an example, when the gradient parameter changes from 2.0 to 5.0, the volume fraction of the ceramic changes from 56.25% to 23.73%, the thermal conductivity increases from 4.45 W/m · K to 6.22 W/m · K, and the temperature of the position point increases from 192.83  to 211.91 K.

Figure 8.

Temperature change of observation point a₃₁ under different gradient parameters.

Figure 9.

The temperature change of observation point a₁₁ under different gradient parameters.

The temperature distribution of the thickness direction with the different gradient parameters is shown in Figure 10 at 2.0 s. It can be seen from the figure that the larger the gradient parameter value is, the higher the temperature is. This is consistent with the influence of gradient parameters on the temperature in x direction. When the observation position in the z direction is close the heat source, the temperature gradually increases. It also can be seen that the closer to the upper surface it is, the greater the rate of temperature change is. This is because the closer to the upper surface it is, the larger the volume fraction of the ceramic is, the lower the thermal conductivity is, and the more dramatic the temperature change. In thickness direction, the volume fraction of ceramics gradually changes, and the heat conduction coefficient also changes gradually. The thermal conductivity has no mutation, so there is no temperature mutation in the thickness direction to avoid the occurrence of heat stress concentration.

Figure 10.

Temperature distribution of the thickness direction with the different gradient parameters (t = 2.0 s).

The temperature varies with the gradient parameter for three observation points on the 0.25 H position from upper surface is shown in Figure 11. From the figure it can be shown that the temperature increases significantly when the gradient parameter ranges from 0 to 3.0, the volume fraction of ceramics changes in the range from 100% to 42.19%, the thermal conductivity changes in the range from 2.09 to 5.22 W/m·K, the thermal conductivity varies in a large range, and the temperature changes significantly. This means the gradient parameters have a significant effect on the temperature distribution in this certain range. When the gradient parameter ranges from 3.0 to 10.0, the volume fraction of ceramics changes in the range from 42.19% to 5.63%, the thermal conductivity changes in the range from 5.22 to 7.20 W/m·K, the thermal conductivity changes in a small range, and the temperature change is relatively small. This means that the gradient parameter has a little effect for the temperature distribution while exceeding a certain range. Explored its causes, the process of heat conduction gradually goes to the pure material during the extreme value of gradient parameters (infinite tends to 0 or ∞). The heat conduction process in the model is changed between pure ceramic heat conduction and pure metal heat conduction with the change of the gradient parameters.

Figure 11.

The temperature varies with the gradient parameter for three observation points on the 0.25 H position from upper surface (t = 2.0 s).

The temperature varies with the gradient parameter for three observation points on the 0.25 H position from lower surface is shown in Figure 12. It can be seen that the peak temperature of the three observation points is below 25 K, and the heat source has a very little effect to the 0.25 H position from the lower surface, this shows the FGM has a good thermal insulation performance. When 3.0 < p <10.0, there is no significant change in temperature. When p = 3.0, the volume fraction of ceramic is about 26.3%, and the volume fraction of metal is about 73.7%. Under this parameter, the thermal insulation property of ceramics and the high strength toughness of metals are well utilized. Therefore, when p = 3.0, the thermal properties and physical properties of FGM will reach the optimal value for the linear heat source.

Figure 12.

The temperature varies with the gradient parameter for three observation points on the 0.25 H position from lower surface (t = 2.0 s).

Conclusion

In this paper, the heat conduction theory of FGM is derived based on the hybrid numerical method, and the temperature distribution in the space-time domain of FGM is solved, the following conclusions are obtained:

1. The change trend of the FGM plate temperature is consistent with the heat source temperature function, showing a linear change. The influence of the heat source on FGM plate is local.

2. The volume fraction of material is determined by gradient parameters in model. The larger the value of the gradient parameter, the higher the peak temperature of the observation point, and the temperature distribution in the model has a positive correlation with the value of the gradient parameter.

3. The gradient parameters have a significant effect for the temperature distribution within a certain range. The heat conduction process in the model is between pure ceramic heat conduction and pure metal heat conduction with the change of the gradient parameters.

4. The value of the gradient parameter determines the physical properties of the material in the model. When the gradient parameter value is 3.0, the thermal properties and physical properties for the linear heat source will reach the optimal value.