Transient Thermoelastic Analysis of Rectangular Plates with Time-Dependent Convection and Radiation Boundaries

: Approximate analytical solutions are presented for the transient thermoelastic problem of rectangular plates with time-dependent convection and radiation boundaries. To include the nonlinear radiation boundary, the whole heating process is divided into several time steps, and a linearized approximation is used to simplify the radiation term for each step. The one-dimensional transient temperature along the thickness direction is solved using the technique of the separation of variables. The displacement and stress solutions are obtained by applying the state-space method to the three-dimensional thermoelasticity equations. The accuracy of the present solutions is demonstrated by comparison with the reported results in the open literature and the ﬁnite element solutions. In the numerical examples, two kinds of thermal boundaries, namely, time-independent convection boundaries and time-dependent convection and radiation boundaries, are considered to show the availability of the present solutions.


Introduction
As fundamental structural elements, plates are extensively applied in building engineering [1][2][3][4].The design and manufacture of the plates needs close attention to the transient thermal and mechanical performance of the structures subjected to complex thermal boundaries.In general, the transient temperature distribution should be obtained before the calculation of the thermal stresses and deformations of the plates [5].
Various numerical and analytical methods have been proposed for the transient heat transfer analysis.The numerical methods, such as the finite element method (FEM) [6] and differential quadrature method [7,8], are well-known tools for solving various transient heat transfer problems.Compared with numerical methods, analytical methods can provide a better insight into the physical significance of different parameters affecting the problem; hence; they are always useful in engineering analysis.The analytical methods used for heat transfer analysis are primarily divided into four categories [9]: the integral transform method [5,10,11], Green's function method [12,13], Laplace transform method (LTM) [14][15][16], and separation of variables [17][18][19][20].Based on the LTM, the one-dimensional (1-D) transient heat conduction in laminated beams with convection boundary was studied by Tanigawa et al. [14].The application of the LTM was restricted in the authors' work since only time-independent (T-I) thermal boundaries were considered.The technique of the separation of variables has been confirmed as effective in dealing with time-dependent (T-D) thermal boundaries.Using this method, de Monte [17,18] studied the 1-D transient heat conduction in composite slabs with convection boundary considering the T-D surrounding temperature.However, the radiation boundary was not considered in these analytical solutions [14,17,18].To include the radiation boundary, Miller and Weaver [19] proposed a linearized method for the radiation term.However, the method is valid only when the surrounding temperature is T-I and the ratio of the surface temperature to the surrounding temperature is larger than 0.75.According to the technique of the separation of variables and orthogonal expansion method, Gu et al. [20] carried out an approximate solution for 1-D temperature distribution in a metallic thermal protection system considering temperature-dependent thermal properties.An advanced linearized method based on Miller and Weaver's work [19] was developed for the radiation term, and an approximate procedure was applied for the T-D surrounding temperature.
For the thermoelastic analysis of plates, most works were carried out based on the approximate two-dimensional (2-D) plate theories [21][22][23][24].Among these approximate theories, the classical plate theory (CPT) is the most commonly used one in engineering analysis due to its simplicity and comprehensibility.Based on the CPT, Yang et al. [25] carried out the stress analysis of an infinite plate with a nonuniformly distributed temperature.Ren et al. [26] carried out the buckling analysis of plates with nonuniformly distributed convection and radiation boundaries.The steady-state heat transfer problem and buckling problem were solved by using the differential quadrature method.Prakash et al. [27] performed an isogeometric analysis for the thermoelastic behaviors of functionally graded (FG) sigmoidal porous plates.
Since the CPT neglects shear deformation effects, it cannot give an accurate analysis for thick plates.To account for the shear deformation effects, researchers have developed diverse shear deformation theories, such as the first-order theory [28][29][30][31][32], third-order theory [33][34][35][36][37], sinusoidal theory [38][39][40][41], and hyperbolic theory [42][43][44].These shear deformation theories can give better predictions on the thermoelastic behaviors of plates in comparison with the CPT.However, they are still only approximations in nature because they take into account shear deformation assumptions, as indicated in the exact threedimensional (3-D) thermoelasticity theory [45].Wang et al. [46] derived exact solutions for laminated Liquid crystal elastomer plates under steady thermal load and mechanical load.A neural network was further trained to predict the deformation of the plates.Aljadani and Zenkour [47] presented a unified form of refined thermoelasticity theory to obtain the field quantities of a rotating/non-rotating half-space with and without the effect of the decay parameter.According to the theory of generalized thermoelasticity [48][49][50] and the memory-dependent derivative, Abouelregal et al. studied the thermomechanical behaviors of a functionally graded piezoelectric rod [51] and an infinite medium [52].Using the Moore-Gibson-Thompson equation, Abouelregal and Alesemi [53] analyzed the thermo-viscoelastic behavior of a fiber-reinforced magnetic viscoelastic solid with temperature-dependent properties.
The above literature review indicates that little attention has been paid to the effects of T-D convection and radiation boundaries on the mechanical performance of plates.Motivated by this fact, we aim to develop approximate analytical solutions to investigate the thermal and mechanical performance of rectangular plates with T-D convection and radiation boundaries.To include the nonlinear radiation term, we divide the whole heating process into several time steps, and use a linearized approximation to simplify the radiation term in each step.The displacement and stress solutions are obtained based on the exact 3-D thermoelasticity theory, which renounces any shear deformation assumptions compared with the approximate 2-D plate theories.In the numerical examples, the transient temperature, displacements, and stresses of a plate are presented considering two kinds of thermal boundaries.

Heat Transfer Analysis
Figure 1 shows a simply supported rectangular plate in the Cartesian coordinate (0 ≤ x ≤ a, 0 ≤ y ≤ b, and 0 ≤ z ≤ H).In the initial stress-free state, the temperature in the plate is uniform and denoted by T 0 .We assume that the plate is subjected to convection and radiation boundary conditions from the top and bottom surfaces by external heat sources with T-D temperatures T ∞t (t) and T ∞b (t), respectively.For simplicity, we assume that the temperature only varies along the thickness direction.

Heat Transfer Analysis
Figure 1 shows a simply supported rectangular plate in the Cartesian coordinate ( ).In the initial stress-free state, the temperature in the plate is uniform and denoted by 0 T .We assume that the plate is subjected to convection and radiation boundary conditions from the top and bottom surfaces by external heat sources with T-D temperatures , respectively.For simplicity, we assume that the temperature only varies along the thickness direction.Before the heat transfer analysis and thermoelastic analysis, some assumptions are established as: i.
The plate is made of a homogeneous, isotropic, and temperature-independent material; ii.The heat convection coefficient and surface emissivity corresponding to the top and bottom surfaces are considered as constants; iii.The plate's deformation falls within the scope of linear elasticity and small strains.

Linearization of Radiation Condition
The plate is subjected to T-D convection and radiation boundary conditions, in which the general radiation condition is expressed by [19] ( ) where q is the heat flux, T ∞ is the surrounding temperature, w T is the temperature at the surface of the plate, σ ( Owing to the radiation condition, the heat transfer problem becomes nonlinear and is hard to exactly solve.Hence, a linear approximation of the radiation condition is required.Based on Miller and Weaver's work [19], Equation (1) can be approximated as It is noted that the error induced by the approximation is neglectable when T ∞ is T- I and 0.75 w T T ∞ > [19].However, the conditions of w T and T ∞ are not satisfied in the present problem, which makes Equation ( 2) not suitable for the present problem.To solve it, we divide the whole heating process into several time steps [12], each of which is assumed to start from zero, i.e., [0, ] t t ∈ Δ .The convection and radiation conditions on the surfaces of the plate within a time step can be written as [20] ( ) It should be noted that by taking j = t and j = b in Equation (3), two equations, i.e., the boundary conditions on the top and bottom surfaces, will be generated.λ is the thermal conductivity, h is the heat convection coefficient, and the factor r is given by [20]  Before the heat transfer analysis and thermoelastic analysis, some assumptions are established as: i.
The plate is made of a homogeneous, isotropic, and temperature-independent material; ii.The heat convection coefficient and surface emissivity corresponding to the top and bottom surfaces are considered as constants; iii.The plate's deformation falls within the scope of linear elasticity and small strains.

Linearization of Radiation Condition
The plate is subjected to T-D convection and radiation boundary conditions, in which the general radiation condition is expressed by [19] where q is the heat flux, T ∞ is the surrounding temperature, T w is the temperature at the surface of the plate, σ (= 5.67 × 10 −8 W/(m 2 K)) is the Stephan-Boltzmann constant, and ε is the surface emissivity.
Owing to the radiation condition, the heat transfer problem becomes nonlinear and is hard to exactly solve.Hence, a linear approximation of the radiation condition is required.Based on Miller and Weaver's work [19], Equation (1) can be approximated as It is noted that the error induced by the approximation is neglectable when T ∞ is T-I and T w /T ∞ > 0.75 [19].However, the conditions of T w and T ∞ are not satisfied in the present problem, which makes Equation ( 2) not suitable for the present problem.To solve it, we divide the whole heating process into several time steps [12], each of which is assumed to start from zero, i.e., t ∈ [0, ∆t].The convection and radiation conditions on the surfaces of the plate within a time step can be written as [20] −λ It should be noted that by taking j = t and j = b in Equation (3), two equations, i.e., the boundary conditions on the top and bottom surfaces, will be generated.λ is the thermal conductivity, h is the heat convection coefficient, and the factor r is given by [20] where T ∞j = 0.

Basic Equations and Homogenization of Boundary Conditions
The transient temperature T(z, t) within a time step is governed by: i. 1-D transient heat transfer equation [12]; where κ = λ/ρc is the thermal diffusivity, ρ is the density, and c is the specific heat.
ii. Simplified convection and radiation conditions on the top and bottom surfaces, as shown in Equation (5); iii.Initial temperature; It is found that the present heat transfer problem contains a homogeneous governing equation (Equation ( 6)) and two non-homogeneous boundary conditions (Equation ( 5)).To homogenize the boundary conditions, a new variable θ(z, t) is introduced as follows: where q(z, t) = (H−z) 2 ζ t T ∞t (t) H(ζ b H+2λ) .Applying Equation (8) to Equations ( 5)- (7) gives It can be seen that the homogenization of the boundary conditions (Equation ( 10)) yields a non-homogeneous governing equation (Equation ( 9)).

Solution of Space-Variable
The corresponding homogeneous equation of Equation ( 9) can be obtained by setting the right side of Equation ( 9) as zero, i.e., Applying the technique of the separation of variables to θ(z, t) gives Substituting Equation ( 13) into Equation ( 12) for each term m yields where β m is the separation constant.Using Equation ( 14) gives the following equation: where only the space-variable Φ m (z) is included.The solution of Equation ( 15) is Substituting Equation ( 16) into Equation (10) gives where According to Cramer's rule, the non-zero solutions of Equation ( 17) exist only if the determinant |C| = 0.The positive roots of |C| = 0 give the eigenvalues Applying each eigenvalue β m to Equation ( 17) obtains A m and B m , where A m can be determined arbitrarily.For example, setting

Solution of Time-Variable
After finding the complete solution of the space-variable Φ m (z), we will solve the time-variable Γ m (t) in this section.To this end, the right side of Equation ( 9), κ∂ 2 q(z, t)/∂z 2 and ∂q(z, t)/∂t, and the initial condition Equation (11), f (z), should be rewritten in the Fourier series form in terms of Φ m (z): Based on the orthogonality of Φ m (z) [16], the unknowns I * m (t), V * m (t), and f * m can be obtained as follows: where 13) and (19) into Equation ( 9) and using Equation (15) gives The initial value Γ m (0) must be given to solve Equation (21).Using Equations ( 11) and ( 13) gives Comparing the third of Equation (19) and Equation ( 22) obtains Γ m (0) = f * m , then the solution of Equation ( 21) is Finally, the transient temperature can be obtained using Equations ( 8), (13), and ( 23)

Calculation Procedure of Temperature Solution
The numerical results of the transient temperature can be obtained in four steps as follows: Step 1: Select an appropriate time interval ∆t i (i where T ∞max denotes the maximum surrounding temperature during the whole heating process [9]. Step 2: The initial temperature in the ith time step is equal to the temperature at the end of the (i − 1)th time step, i.e., F i (z) = T i−1 (z, ∆t i−1 ).
Step 3: The transient temperature T i (z, t i ) in the ith time step can be obtained based on the method introduced in Sections 2.1-2.4.

Thermoelasticity Equations
The thermoelastic analysis of the plate is within the framework of 3D thermoelasticity theory.Thus, the deformation and stress states are governed by [45,54]: i. Geometrical relations; iii.Thermoelastic constitutive relations; iv. Simply supported edge conditions; v. Surface stress conditions; In Equations ( 25)-( 29) {σ x σ y σ z τ xy τ xz τ yz } are stress components, {ε x ε y ε z γ xy γ xz γ yz } strain components, u v w displacement components, and where α is the thermal expansion coefficient, E is Young's modulus, and µ is Poisson's ratio.
Applying the state-space method [45] to Equations ( 25)- (27) gives where The stress components σ x , σ y , and τ xy can be solved by

General Displacement and Stress Solutions
Because the plate is simply supported at its four edges, the displacement and stress solutions can be assumed as where a r = rπ/a and b s = sπ/b.
Similarly, the thermoelastic term P(z, t) in Equations ( 31) and ( 33) can be expanded as where P r,s (z, t) = 4 ab a 0 b 0 P(z, t) sin(a r x) sin(b s y)dxdy.Applying Equations ( 34) and (35) to Equations (33) gives According to Equations ( 34) and (36), we can find that the edge conditions (see Equation ( 28)) are satisfied.
Applying Equations ( 34) and (35) to Equation (31) gives where Based on the matrix theory, the solution to Equation ( 37) is [45] P(z, t) = M(z)P(0, t) + N(z, t) where To calculate M(z) and N(z, t), the eigenvalues of A should be taken into consideration.Let κ i (i = 1, 2, . .., 6) be the eigenvalues of A and J i (i = 1, 2, . .., 6) be the corresponding eigenvectors, then there must exist a matrix K = [J 1 J 2 • • • J 6 ] and its inverse matrix K −1 which can transform A into a diagonal matrix.In this case, we have Therefore

Determination of Unknowns
By applying Equations ( 34) and (35) to Equation ( 29), the surface stress conditions can be rewritten as Taking z = H into Equation (39) yields After eliminating U r,s (H, t), V r,s (H, t), and W r,s (H, t) from Equation ( 44), we have where M ij is the element of the matrix M(H) and N i is the element of the vector N(H, t).
Finally, the displacement and stress solutions are obtained from Equations ( 34) and (36).

Numerical Results and Discussion
In the following numerical calculation, the first 40 eigenvalues (i.e., β 1 , β 2 , . .., β 40 ) are used to obtain the temperature solution, and the first 25 × 25 terms (i.e., r = s = 1, 2, . .., 25) are used to obtain the displacement and stress solutions.Then, the comparison of the thermal bending behavior is carried out.Fogang [56] provided the thermal bending analysis of a rectangular plate under linear temperature distribution across the thickness.Here, we use the present method to estimate the central deflection of the plate, and compare our results with the published ones [56].The deflection is expressed in the following nondimensional form:

(
) where t T Δ and b T Δ are the temperature increments at the top and bottom surfaces of the plate, respectively.Table 1 shows the comparative results for different aspect ratios a/b and different Poisson's ratios.The present results match well with those reported by Fogang [56].The maximum error is less than 0.5%.Then, the comparison of the thermal bending behavior is carried out.Fogang [56] provided the thermal bending analysis of a rectangular plate under linear temperature distribution across the thickness.Here, we use the present method to estimate the central deflection of the plate, and compare our results with the published ones [56].The deflection is expressed in the following nondimensional form: where ∆T t and ∆T b are the temperature increments at the top and bottom surfaces of the plate, respectively.Table 1 shows the comparative results for different aspect ratios a/b and different Poisson's ratios.The present results match well with those reported by Fogang [56].The maximum error is less than 0.5%.

Purely Convection Boundary and T-I Surrounding Temperature
Consider a rectangular steel plate with the sizes given by the following: a = 0.8 m, b = 0.6 m, and H = 0.1 m.The material properties of the plate are as follows [45]: λ = 54.33W/(m • K), ρ = 7850 kg/m 3 , c = 440 J/(kg • K), ε = 0.7, E = 210 GPa, µ = 0.3, and α = 1.216 × 10 −5 K −1 .The initial temperature is fixed at T 0 = 300 K.The heat convection coefficients of the top and bottom surfaces are h t = 800 W/(m 2 • K) and h b = 600 W/(m 2 • K), respectively.The bottom surrounding temperature is T-I and fixed at T ∞b = 300 K, while the top surrounding temperature T ∞t is variable.
Consider that the plate is subjected to purely convection boundary conditions.The top surrounding temperature is T-I and fixed at T ∞t = 375 K. Owing to the omission of the radiation boundary, the whole heating process need not be divided into several time steps.A 3-D FEM analysis is carried out for comparison by use of the software ABAQUS.The transient heat transfer and thermoelastic behaviors are simulated by using the 8-node element DC3D8 and C3D8R, respectively, and by meshing the plate with 80, 60, and 20 elements in the x-, y-, and z-directions, respectively.
Figure 3 shows the comparative results among the time history of the temperature obtained from the present method, the FEM, and the LTM [11] at z = 0, z = 0.05 m, and z = 0.1 m.A good agreement can be found among the three methods.The temperature changes quickly in the early stage and tends to be stable after about t = 900 s.It is also noticed that the higher the position is, the earlier the temperature starts to rise.This is because the heat flows from the top surface (z = 0.1 m) to the bottom surface (z = 0).
Figure 4 shows the comparative results among the temperature distributions obtained from the present method, the FEM, and the LTM [11] at t = 10 s, t = 100 s, and t = 1000 s.It is noticed that the accuracy of the FEM solution improves with time, since the FEM simulation cannot provide a highly accurate prediction when the temperature changes quickly at the early stage.It is also seen that the temperature distribution progressively changes from nonlinearity to linearity with time.This is because of the homogeneity and temperature independency of the material properties, resulting in linear heat conduction in the equilibrium stage.
tained from the present method, the FEM, and the LTM [11] at t = 10 s, t = 100 s, and t = 1000 s.It is noticed that the accuracy of the FEM solution improves with time, since the FEM simulation cannot provide a highly accurate prediction when the temperature changes quickly at the early stage.It is also seen that the temperature distribution progressively changes from nonlinearity to linearity with time.This is because of the homogeneity and temperature independency of the material properties, resulting in linear heat conduction in the equilibrium stage.tained from the present method, the FEM, and the LTM [11] at t = 10 s, t = 100 s, and t = 1000 s.It is noticed that the accuracy of the FEM solution improves with time, since the FEM simulation cannot provide a highly accurate prediction when the temperature changes quickly at the early stage.It is also seen that the temperature distribution progressively changes from nonlinearity to linearity with time.This is because of the homogeneity and temperature independency of the material properties, resulting in linear heat conduction in the equilibrium stage.Figure 5 shows the through-thickness distributions of u, w, and σ x at x = 0.2 m, y = 0.15 m for different times t = 10 s, t = 100 s, and t = 1000 s.It is seen that the displacement and stress distributions obtained from the present method and the FEM have a good agreement.Figure 5a-c show that the deformations and stresses are increased with time because the thermal load is increased with time.The displacements u and w vary slightly with the coordinate z at t = 10 s compared with those at t = 1000 s. Figure 5c shows that the distribution of σ x changes from nonlinearity to linearity with time, which is similar to the variation of the temperature distribution.

Convection and Radiation Boundaries and T-D Surrounding Temperature
Consider that the plate is subjected to convection and radiation boundary conditions.The top surrounding temperature is T-D and expressed as where the coefficient η represents the rate of temperature variation.Figure 6 shows the time history of T ∞t for different coefficients η = 0.005, η = 0.01, and η = 0.05.In the following analysis, the interval of each time step is uniform and is taken as ∆t i = 2 s.
0.15 m for different times t = 10 s, t = 100 s, and t = 1000 s.It is seen that the displacement and stress distributions obtained from the present method and the FEM have a good agreement.Figure 5a-c show that the deformations and stresses are increased with time because the thermal load is increased with time.The displacements u and w vary slightly with the coordinate z at t = 10 s compared with those at t = 1000 s. Figure 5c shows that the distribution of x σ changes from nonlinearity to linearity with time, which is similar to the variation of the temperature distribution.

Convection and Radiation Boundaries and T-D Surrounding Temperature
Consider that the plate is subjected to convection and radiation boundary conditions.The top surrounding temperature is T-D and expressed as where the coefficient η represents the rate of temperature variation.Figure 6 shows the time history of t T ∞ for different coefficients η = 0.005, η = 0.01, and η = 0.05.In the follow- ing analysis, the interval of each time step is uniform and is taken as    Figure 7 shows the time history of temperature at z = 0 and z = 0.1 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It can be seen that the present results using ∆t i = 2 s match well with the FEM ones.The temperature rises faster with η, since the top surrounding temperature T ∞t rises faster with η (see Figure 6).Figure 7 shows the time history of temperature at z = 0 and z = 0.1 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It can be seen that the present results using 2 s i t Δ = match well with the FEM ones.The temperature rises faster with η, since the top surrounding temperature t T ∞ rises faster with η (see Figure 6).Figure 8 shows the temperature distributions at t = 100 s and t = 1000 s for different coefficients η = 0.005, η = 0.01, and η = 0.05.It is seen from Figure 8a,b that the temperature distributions are nonlinear at t = 100 s, but are nearly linear at t = 1000 s, which indicates that the temperature in the plate tends to be stable at about t = 1000 s.It is also noticed that, at t = 100 s, the FEM solution is less accurate with the increase in η; however, a similar phenomenon does not occur at t = 1000 s.This is because the temperature in the plate rises faster with η at the early stage, and tends to be stable at about t = 1000 s. Figure 8 shows the temperature distributions at t = 100 s and t = 1000 s for different coefficients η = 0.005, η = 0.01, and η = 0.05.It is seen from Figure 8a,b that the temperature distributions are nonlinear at t = 100 s, but are nearly linear at t = 1000 s, which indicates that the temperature in the plate tends to be stable at about t = 1000 s.It is also noticed that, at t = 100 s, the FEM solution is less accurate with the increase in η; however, a similar phenomenon does not occur at t = s.This is because the temperature in the plate rises faster with η at the early stage, and tends to be stable at about t = 1000 s.    Figure 9 shows the through-thickness distributions of w and τ xy at x = 0.48 m, y = 0.48 m, and t = 100 s for different coefficients η = 0.005, η = 0.01, and η = 0.05. Figure 9a,b indicate that the maximum absolute values of w and τ xy are increased with η.The slopes of the curves of w and τ xy are both increased with η since the temperature gradient is increased with η (see Figure 8a).It is also noticed from Figure 9a that w is variable with the coordinate z.This demonstrates that the approximate 2-D plate theories, which assume that w is invariable through the thickness direction, may generate errors when estimating the thermomechanical responses of the plate.It is also found that w increases quickly and tends to be stable after about t = 900 s; however, xy τ still has an obvious increase at t = 1000 s. Figure 10a indicates that the long-term deformations (i.e., deformations at t→∞) of the plate do not vary with η.To show the effects of radiation heat transfer on the mechanical performance of the plate, two kinds of thermal boundaries are considered: (i) purely convection boundary and (ii) convection and radiation boundaries.Figure 11 shows the through-thickness dis- Figure 10 shows the time history of w and τ xy at x = 0.48 m, y = 0.48 m, and z = 0.05 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It is seen that w and τ xy increase quickly with η.It is also found that w increases quickly and tends to be stable after about t = 900 s; however, τ xy still has an obvious increase at t = 1000 s. Figure 10a indicates that the long-term deformations (i.e., deformations at t→∞) of the plate do not vary with η.To show the effects of radiation heat transfer on the mechanical performance of the plate, two kinds of thermal boundaries are considered: (i) purely convection boundary and (ii) convection and radiation boundaries.Figure 11  To show the effects of radiation heat transfer on the mechanical performance of the plate, two kinds of thermal boundaries are considered: (i) purely convection boundary and (ii) convection and radiation boundaries.Figure 11 shows the through-thickness distributions of w and τ xy at x = 0.48 m, y = 0.48 m, and t = 1000 s for η = 0.01.It is noted that to reduce convection effects, the heat convection coefficients of the top and bottom surfaces are taken as h t = h b = 10 W/(m 2 • K).It is seen that distributions of w and τ xy are greatly affected by radiation heat transfer.The values of w and τ xy with radiation effects are larger than those without radiation effects because radiation heat transfer will generate a larger thermal load.

Conclusions
This work presents approximate analytical solutions for the transient thermoelastic problem of rectangular plates with T-D convection and radiation boundaries.By using a linearized approximation to simplify the radiation term, 1-D temperature distribution is obtained.The displacement and stress distributions are obtained based on the exact 3-D thermoelasticity theory.In the numerical examples, the temperature, displacement, and stress fields in a rectangular plate are studied considering two kinds of thermal boundaries.It is concluded that: i.
The comparison with the reported results and the FEM results indicates that the present method can give a satisfactory prediction for the thermoelastic responses in the plate.ii.The temperature, displacements, and stresses increase quickly in the early stage but tend to eventually stabilize.iii.The temperature distribution progressively changes from nonlinearity to linearity with time.iv.The temperature, displacements, and stresses increase quickly with the increase in the temperature coefficient η.However, the long-term deformations of the plate do not vary with η. v.The radiation heat transfer has a significant effect on the displacement and stress distributions.

Conclusions
This work presents approximate analytical solutions for the transient thermoelastic problem of rectangular plates with T-D convection and radiation boundaries.By using a linearized approximation to simplify the radiation term, 1-D temperature distribution is obtained.The displacement and stress distributions are obtained based on the exact 3-D thermoelasticity theory.In the numerical examples, the temperature, displacement, and stress fields in a rectangular plate are studied considering two kinds of thermal boundaries.It is concluded that: i.
The comparison with the reported results and the FEM results indicates that the present method can give a satisfactory prediction for the thermoelastic responses in the plate.ii.The temperature, displacements, and stresses increase quickly in the early stage but tend to eventually stabilize.iii.The temperature distribution progressively changes from nonlinearity to linearity with time.iv.The temperature, displacements, and stresses increase quickly with the increase in the temperature coefficient η.However, the long-term deformations of the plate do not vary with η. v.
The radiation heat transfer has a significant effect on the displacement and stress distributions.

Figure 1 .
Figure 1.A rectangular plate with T-D convection and radiation boundaries.
Stephan-Boltzmann constant, and ε is the surface emissivity.

Figure 1 .
Figure 1.A rectangular plate with T-D convection and radiation boundaries.

4. 1 .
Validation Study As no published results for the rectangular plates under current consideration are available in the open literature, two examples concerning the transient heat transfer and thermal bending behaviors are used to validate the present methods by comparing our results with the published ones.First, the comparison of the transient heat transfer behavior is carried out.Monds and McDonald [55] built two mathematical models, namely, the finite-length scale model and semi-infinite model, to estimate the 1-D transient temperature in human skin during a simulated fire.The accuracy of the models was checked by comparison with the experimental results.In their work, the material properties and thickness of the skin simulant are the following: λ = 0.97 W/(m • • C), ρ = 1877 kg/m 3 , c = 1205 J/(kg • • C), and H = 19.66mm.The temperatures at the top and bottom surfaces of the skin simulant are the following: T t = 0.0013t 2 − 0.0696t + 21.757 • C and T b = −0.0093t 2 + 2.1531t + 60.255 • C. Here, we use the present method to estimate the temperature variations at z = 11.12 mm and compare our results with the published ones [55].The comparative results are shown in Figure 2. The present results match well with those of the finite-length scale model, but have a slight difference from the experimental ones when t > 100 s.The results of the semi-infinite model are more inaccurate compared with the present method and the finite-length scale model.Buildings 2023, 13, x FOR PEER REVIEW 10 of 19

Figure 3 .
Figure 3.Time history of temperature at z = 0, z = 0.05 m, and z = 0.1 m considering convection boundary and T-I surrounding temperature.

Figure 3 .
Figure 3.Time history of temperature at z = 0, z = 0.05 m, and z = 0.1 m considering convection boundary and T-I surrounding temperature.

Figure 3 .
Figure 3.Time history of temperature at z = 0, z = 0.05 m, and z = 0.1 m considering convection boundary and T-I surrounding temperature.

Figure 4 .
Figure 4. Through-thickness temperature distributions at t = 10 s, t = 100 s, and t = 1000 s considering convection boundary and T-I surrounding temperature.

Figure 5 .
Figure 5. Through-thickness displacement and stress distributions at x = 0.2 m, y = 0.15 m considering convection boundary with T-I surrounding temperature: (a) displacement u; (b) displacement w; (c) stress x σ .

Figure 7
Figure7shows the time history of temperature at z = 0 and z = 0.1 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It can be seen that the present results using2 s i t Δ =

Figure 7 .
Figure 7. Time history of temperature for different coefficients η = 0.005, η = 0.01, and η = 0.05 considering convection and radiation boundaries and T-D surrounding temperature: (a) temperature T at z = 0; (b) temperature T at z = 0.1 m.

Figure 7 .
Figure 7. Time history of temperature for different coefficients η = 0.005, η = 0.01, and η = 0.05 considering convection and radiation boundaries and T-D surrounding temperature: (a) temperature T at z = 0; (b) temperature T at z = 0.1 m.

Figure 9
Figure 9 shows the through-thickness distributions of w and xy τ at x = 0.48 m, y = 0.48 m, and t = 100 s for different coefficients η = 0.005, η = 0.01, and η = 0.05. Figure 9a,b indicate that the maximum absolute values of w and xy τ are increased with η.The slopes

Figure 10 shows
Figure10shows the time history of w and xy τ at x = 0.48 m, y = 0.48 m, and z = 0.05 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It is seen that w and xy τ increase quickly with η.It is also found that w increases quickly and tends to be stable after about t = 900 s; however, xy

Figure 10 shows
Figure10shows the time history of w and xy τ at x = 0.48 m, y = 0.48 m, and z = 0.05 m for different coefficients η = 0.005, η = 0.01, and η = 0.05.It is seen that w and xy τ increase quickly with η.It is also found that w increases quickly and tends to be stable after about t = 900 s; however, xy τ still has an obvious increase at t = 1000 s.Figure10aindicates that the long-term deformations (i.e., deformations at t→∞) of the plate do not vary with η.

Figure 11 .
Figure 11.Effects of radiation heat transfer on the through-thickness displacement and stress distributions at x = 0.48 m, y = 0.48 m, and t = 1000 s for η = 0.01: (a) displacement w; (b) stress xy τ .

Funding:Figure 11 .
Figure 11.Effects of radiation heat transfer on the through-thickness displacement and stress distributions at x = 0.48 m, y = 0.48 m, and t = 1000 s for η = 0.01: (a) displacement w; (b) stress τ xy .

Table 1 .
Central deflection for different aspect ratios a/b and different Poisson's ratios.

Table 1 .
Central deflection for different aspect ratios a/b and different Poisson's ratios.
Author Contributions: Conceptualization, Z.Z., Z.X. and X.S.; methodology, Z.Z. and Y.S.; software, Z.Z.; validation, Z.Z., Y.S. and W.Q.; formal analysis, Z.X. and W.Q.; investigation, Z.Z. and Y.S.; resources, Z.X., W.Q. and X.S.; writing-original draft preparation, Z.Z.; writing-review and editing, Z.Z. and Z.X.; project administration, Z.Z. and Z.X.; funding acquisition, W.Q. and X.S.All authors have read and agreed to the published version of the manuscript.Funding: This work is financially supported by the National Natural Science Foundation of China (Grant No. 52208395) and the Nantong Science and Technology Plan Project (MS2021030).