Strong and Weak Formulations of a Mixed Higher-Order Shear Deformation Theory for the Static Analysis of Functionally Graded Beams under Thermo-Mechanical Loads

: The strong and weak formulations of a mixed layer-wise (LW) higher-order shear deformation theory (HSDT) are developed for the static analysis of functionally graded (FG) beams under various boundary conditions subjected to thermo-mechanical loads. The material properties of the FG beam are assumed to obey a power-law distribution of the volume fractions of the constituents through the thickness of the FG beam, for which the e ﬀ ective material properties are estimated using the rule of mixtures, or it is directly assumed that the e ﬀ ective material properties of the FG beam obey an exponential function distribution along the thickness direction of the FG beam. The results shown in the numerical examples indicate that the mixed LW HSDT solutions for elastic and thermal ﬁeld variables are in excellent agreement with the accurate solutions available in the literature. A parametric study related to various e ﬀ ects on the coupled thermo-mechanical behavior of FG beams is carried out, including the aspect ratio, the material-property gradient index, and di ﬀ erent boundary conditions. approach 3, they are directly obtained using the Lagrange multipliers. Implementation of the three approaches shows that the accuracy of these


Introduction
Functionally graded (FG) structures are emerging composite structures, for which material properties can be designed to gradually and smoothly vary through their physical domains. FG structures can be formed by mixing two-phase materials with pre-designed single-and bi-directional distributions of the volume fractions of the constituents through the thickness coordinate and the axial-thickness surface of the FG structures, respectively [1][2][3][4][5]. Due to their graded material properties, FG structures can prevent delamination and stress concentration phenomena, usually occurring at the interfaces between adjacent layers in the cases of laminated fiber-reinforced composite (FRC) structures as a result of the material properties suddenly changing at these locations. Some FG thermoelectric devices were developed to enhance their thermal stability resistance. For example, FG piezoelectric materials were developed for broad-band ultrasonic transducers [6] and FG composite electrodes were developed for solid oxide fuel cells [7]. The functional gradation principle was used to develop an artificial biomaterial for knee joint replacement [8]. Some FG materials were applied to optimize the thermal, wear, and corrosion properties of metallic and ceramic materials [9]. For example, FG thermal barrier coating materials were deposited onto Cu substrates to improve the interface fracture toughness between the coatings and the substrates. FG TiC-TiN thin films were used to improve the wear resistance of cutting tool alloys with good adhesion to the substrate material. FG materials were obtained using the DQ method, respectively. These quasi-3D solutions can provide a standard by which to assess the accuracy and convergence rates of assorted advanced and refined ESLTs and LW theories available in the literature.
After a close literature survey, it was found that most of the LW theories were applied to the mechanical analyses of laminated FRC beams and were rarely applied to the mechanical analyses of FG beams. Due to the excellent performance of the mixed LW HSDT [22,23], as mentioned above, the authors extend it to the current coupled thermo-mechanical analysis of FG beams with various boundary conditions. In the formulation, the displacement components are expanded as an LW third-order polynomial function through the thickness direction of the beam, and the displacement continuity conditions at the interfaces between adjacent layers are introduced in the potential energy functional using the Lagrange multiplier method, such that they are satisfied in a variational form. The Lagrange multipliers are the exact transverse stress components. Based on the stationary principle of the extended potential functional, the strong and weak formulations of the mixed LW HSDT are derived. The Navier-type analytical solutions based on the strong formulation and the finite element (FE) solutions based on the weak formulation for the static analysis of FG beams subjected to thermo-mechanical loads are obtained. The accuracy and convergence rates of these analytical and finite element solutions are validated by comparing them with the quasi-3D solutions available in the literature. Some effects on the coupled thermo-mechanical behavior of FG beams are conducted, including the aspect ratio, the material-property gradient index, and different boundary conditions.

Effective Material Properties
In this work, the rule of mixtures [32] is used to estimate the effective material properties of the FG beam and is described as follows: According to the rule of mixtures, the through-thickness distributions of the effective material properties of the FG beam can be written in the following form, where Γ c and Γ m represent the volume fractions of the ceramic and metal materials of the constituents of the FG beam, respectively, such that Γ m + Γ c = 1. F can be one of the engineering constants, including Young's modulus E, Poisson's ratio υ, thermal expansion coefficient α, and mass density ρ.
The subscripts m and c represent a metal material and a ceramic material, respectively.

The Strong Formulation and Its Application
In this section, the authors develop the strong formulation of a mixed LW HSDT for a static analysis of FG beams under various boundary conditions subjected to thermo-mechanical loads, where the material properties of the FG beam are considered to be thickness dependent. The configuration and coordinates of the FG beam are shown in Figure 1, where h and L represent the thickness and the length of the FG beam, respectively. In the analysis, the FG beam is artificially divided into n l layers, and the thickness of each individual layer constituting the beam is h m (m = 1 − n l ), such that n l m=1 h m = h.
In the mixed LW HSDT, the displacement field for a typical individual layer is given as follows: w (m) (x, z m ) = w where m = 1, 2, . . . , n l , and u  In the mixed LW HSDT, the displacement field for a typical individual layer is given as follows: where m = 1, 2, …, nl, and ( ) The displacement components in the y direction are taken to be zero.
According to the perfect bonding assumptions at the interfaces between adjacent layers, the corresponding displacement continuity conditions at these locations are given as where k = 1, 2,…, (nl−1). The strain-displacement relationship is given as  According to the perfect bonding assumptions at the interfaces between adjacent layers, the corresponding displacement continuity conditions at these locations are given as where k = 1, 2, . . . , (n l − 1). The strain-displacement relationship is given as γ (m) where the commas denote the derivative of the suffix variable, and the remaining strains are zeroes, xy . The stress-strain relationship for an orthotropic material in thermal environment is given as where ∆T denotes the temperature change, measured at a room temperature of 300 K. Q (m) 3 in which the subscripts 1, 2, and 3 denote the principle axes of the material properties, and E, α and υ and represent the Young's modulus, the thermal expansion coefficients, and Poisson's ratio, respectively. For isotropic materials, these stiffness coefficients will be reduced to Q (m) Note that these engineering constants E, α and υ in the analysis are considered to be dependent from the thickness of the FG beam, where ∆T is a function of x and z, i.e., ∆T(x, z).
The governing equations and associated boundary conditions are derived using the stationary principle of minimum potential energy combined with the Lagrange multiplier method, in which the displacement continuity conditions at the interfaces between adjacent layers given in Equations (4) and (5) are multiplied by the Lagrange multipliers and are then substituted into the potential energy functional as the constraints, such that the extended potential energy functional of the n l -layered FG beam can be given as follows: where σ (m) x and τ z ) induced at the interfaces between adjacent layers, respectively. q(x) is the external load applied on the top surface of the FG beam, where its positive direction is defined to be downward.
Equations (31)-(40), associated with a set of boundary conditions (i.e., Equations (23)-(30)), can be composed of a well-post boundary value problem, for which the Navier-type analytical solutions of the elastic and thermal variables induced in simply supported, multi-layered FG beams can be obtained using the Fourier series expansion method, for which the detailed solution process is given in Appendix B.

The Weak Formulation and Its Application
To develop the weak formulation of the mixed LW HSDT, the authors perform the first-order variational operation with respect to the potential energy functional of the multi-layered FG beams, as follows: where n e denotes the number of nodes used for a typical element.
The primary variables are expressed as where n d denotes the number of nodes, and n d = 2, 3, and 4 is used for the linear, quadratic, and cubic elements, respectively. The strain and stress components, the Lagrange multipliers, and the displacement continuity conditions are expressed in the matrix form as follows: where Substituting Equations (42) and (43) into Equation (41) yields where The finite element solutions of the nodal displacement components and the nodal Lagrange multipliers can be obtained by solving Equation (44). The secondary field variables can thus be obtained using the determined primary nodal field variables.

Boundary Conditions
There are four boundary conditions considered in the following numerical examples, i.e., C-C, C-S, S-S, and C-F boundary conditions, which are given as follows: For the C-C boundary conditions, For the C-S boundary conditions, For the S-S boundary conditions, For the C-F boundary conditions,

Numerical Examples
Based on the above-mentioned strong and weak formulations, the authors can obtain the Navier-type analytical solutions for static problems of simply supported, multi-layered FG beams subjected to thermo-mechanical loads, as well as the FE numerical solutions for those of multi-layered FG beams under various boundary conditions, respectively. Some numerical examples are examined and discussed in the following sections.

Mechanical Loads
In this work, the static behavior of an FG beam under various boundary conditions subjected to a uniformly distributed mechanical load (i.e., q(x) = q 0 ) applied on the top surface of an FG beam is considered. The FG beam is composed of a two-phase composite material by mixing a metal material (aluminum) and a ceramic material (alumina) according to a power-law distribution of the volume fractions of the constituents through the thickness direction of the FG beam. The material properties of the aluminum and alumina materials are E m = 70 GPa, υ m = 0.3, and E c = 380 GPa, υ c = 0.3, respectively, where the subscripts m and c represent the metal (aluminum) and ceramic (alumina) materials. For comparison purposes, a set of dimensionless variables is defined as the same forms as those used in Thai and Vo [16], and can be given as follows: The applied external load is expanded as the single Fourier series in the x direction as follows: where qm = 4q 0 / m π , whenm is an odd integer, and qm = 0, whenm is an even integer. In this example,m is taken asm = 1, 3, 5, . . . , 39 for a uniform mechanical load. The volume fraction of the ceramic material Γ c (z) is defined as Γ c (z) = [(1/2) + (z/h)] κ p , where κ p denotes the material-property gradient index of a power-law-type material model. The effective material properties are estimated using the rule of mixtures, for which a formula is given as Equation (1).
In the mixed LW HSDT, the transverse shear and normal stress components induced in the loaded FG beams can be obtained using three different approaches. In approach 1, they are obtained using constitutive equations when the displacement components are determined. In approach 2, they are obtained using stress equilibrium equations when the displacement components are determined. In approach 3, they are directly obtained using the Lagrange multipliers. Table 1 shows the convergent study for analytical solutions of the stress and displacement components induced in a simply supported, FG beam subjected to a uniform mechanical load, where L/h = 5 and 20, and L = 5 m; κ p = 5; n l = 4, 8, and 16 for the LW1, and n l = 2, 4, and 8 for the LW2 and LW3 theories. LW1, LW2, and LW3 represent the mixed LW FSDT, the mixed LW second-order shear deformation, and the mixed LW TSDT, respectively. It can be seen in Table 1 that the convergent solutions are obtained when n l is taken to be 16, 8, and 4, for the LW1, LW2, and LW3 theories, respectively. The convergence rates of LW1, LW2, and LW3 are in the following order: LW3 > LW2 > LW1, where the symbol ">" indicates a fast convergence rate. The convergent results for the LW1, LW2, and LW3 theories are also compared with the solutions obtained using the HSDT [33], RSDT [34], sinusoidal shear deformation theory (SSDT) [35], TSDT [36], exponential shear deformation theory (ESDT) [37], and classical beam theory (CBT) [16]. The results show that the convergent solutions of the displacement and in-plane stress components obtained using the mixed LW HSDT are in excellent agreement with those solutions obtained using the HSDT, RHSDT, SSDT, TSDT, and ESDT. The CBT fails to provide satisfactory results for the analysis of FG beams. It is well known that the transverse shear and normal stress components obtained using the stress equilibrium equation approach (i.e., approach 2) are significantly more accurate than those obtained using the constitutive equation approach (i.e., approach 1), although the solution process for the former is more complicated and more time consuming than that of the latter. In the LW HSDTs, the transverse shear stress components can also be directly obtained using the Lagrange multipliers (i.e., approach 3), which are the primary variables, without any additional processing procedures, and these components obtained using approach 3 closely agree with those obtained using approach 2. In most of the beam theories available in the literature, the transverse shear and normal stress components were obtained using the constitutive equation approach (i.e., approach 1), such that these components predicted by various theories are different from each other. On the basis of the convergent solutions for LW3, the relative errors for the transverse shear stress components obtained using the HSDT, RHSDT, SSDT, TSDT, and ESDT are 9.42%, 7.62%, 3.71%, 7.96%, and 0.19%, respectively, for thick beams (L/h = 5). The accuracy of the transverse shear stress components for various theories is in the following order: ESDT > SSDT > (RHSDT, TSDT) > HSDT, where the symbol ">" represents greater accuracy. In order to expand the applicable range of the mixed LW HSDT to the analysis of FG beams under different boundary conditions, rather than under simply supported boundary conditions only, a weak formulation-based finite element method is developed in this work. As the results in Table 1 show that the performance of LW3 is superior to that of LW2 and LW1, an LW3-based quadratic finite element method is developed for the thermo-mechanical analysis of FG beams under combinations of simply supported, free, and clamped boundary conditions. Table 2 shows a convergent study for the LW3-based quadratic FEM solutions of stress and displacement components induced in the FG beam under various boundary conditions, including the C-C, C-S, S-S, and C-F boundary conditions, subjected to a uniform mechanical load. It can be seen in Table 2 that the convergent solutions for the stress and displacement components induced in the FG beam can be obtained using an element mesh (n x × n z ) = (32 × 16) for a wide range of L/h > 5, where n x and n z denote the numbers of elements used in the x and z directions, respectively. For a moderately thick FG beam (L/h = 10) and on the basis of an allowable relative error of 1%, the convergent solutions of displacement, in-plane stress, and transverse stress components can be obtained when (n x × n z ) = (8 × 16), (16 × 16), and (32 × 16) element meshes are used, respectively. The values of the displacement components u x and u z at some specific positions for assorted boundary conditions are in the following order: C-C < C-S < S-S < C-F, where the symbol "<" indicates a less displacement and greater gross beam stiffness. Figures 2 and 3 show variations in the through-thickness distributions of various stress and displacement components induced in a simply supported FG beam under a uniform mechanical load with different L/h ratios and different values of material-property gradient indices, respectively, where the LW3 theory with n l = 32 is used. In Figure 2, the length-to-thickness ratio is taken as L/h = 5, 10, and 20, and the material-property gradient index is taken as κ p = 3. In Figure 3, they are L/h = 10, and κ p = 0, 0.2, 5, and ∞. where the LW3 theory with nl = 32 is used. In Figure 2, the length-to-thickness ratio is taken as L/h = 5, 10, and 20, and the material-property gradient index is taken as p κ = 3. In Figure 3, they are L/h = 10, and p κ = 0, 0.2, 5, and ∞.       A set of dimensionless variables is defined as follows: It can be seen in Figures 2b and 3b that the in-plane normal stress component σ x induced in the FG beam increases when the FG beam becomes thinner, where κ p = 3. It can also be seen that the through-thickness distribution of σ x induced in the FG beam appears to be a higher-order polynomial function, κ p = 3 and 5 rather than a linear function induced in an isotropic homogeneous beam (i.e., κ p = 0). The results in Figure 2c and Figure 3c show that the through-thickness distribution of transverse shear stress component τ xz appears to be a parabolic function for the isotropic homogeneous beam (i.e., κ p = 0), where its peak value occurs in the mid-plane, and it appears to be a higher-order polynomial function with a peak value occurring at a position above half the beam height for the FG beam (i.e., κ p 0. Referring to the relationships of the dimensionless displacement component w to its dimensional counterpart w, it can be seen in Figures 2a and 3a that the displacement component w increases when the value of κ p becomes greater, which indicates that the volume fraction of the metal material becomes greater than that of the ceramic material. In Figure 2c,d and Figure 3c,d, it is shown that the traction conditions on the top surface of the FG beam are exactly satisfied.

Thermal loads
In this section, the authors investigate the stress and deformation components induced in an FG beam under combinations of simply supported, free, and clamped boundary conditions subjected to a sinusoidally distributed thermal load using the strong and weak formulations of the mixed LW HSDTs. For comparison purposes, a numerical example used in Lü et al. [31] is adopted to validate the accuracy and convergence rate of the strong formulation-based analytical solutions and the weak formulation-based finite element solutions. In this example, the length and thickness of the FG beam are taken as L = 0.1 m and h = 0.01 m, respectively. The Young's modulus, Poisson's ratio, thermomechanical coefficients, and temperature changes are defined as where E 0 ,, Q α0 , and T 0 denote a reference modulus, a reference thermomechanical coefficient, and a reference temperature change, respectively; Q α0 = 10 −6 E 0 /(1 − υ), and T 0 = 100 K κ e stands for the material-property gradient index of an exponential function-type material model. A set of dimensionless displacement and stress components is defined as having the same form as that used in Lü et al. [31] and can be given as follows: Table 3 shows the convergent study for analytical solutions of thermal stress and thermal displacement components induced in a simply supported, FG beam subjected to a sinusoidally distributed thermal load, where κ e = 230 and 100; h/L = 0.1 and 0.2, and L = 0.1 m; n l = 16, 32, and 64 for the LW1 theory, n l = 8, 16, and 32 for the LW2 theory, and n l = 4, 8, and 16 for the LW3 theory. It can be seen in Table 3 that the convergent solutions are obtained when n l is taken to be 32, 16, and 8, for the LW1, LW2, and LW3 theories, respectively. Again, the convergence rates of the LW1, LW2, and LW3 are in the following order: LW3 > LW2 > LW1, where the symbol ">" indicates a fast convergence rate. The convergent results of the LW1, LW2, and LW3 theories are also compared with the solutions obtained using a 13-node state space-based differential quadrature method (SSDQM) [31] and a bilinear rectangular plane element method (BRPEM) with an element mesh (n x × n z ) = (200 × 20) [31]. The results show that the convergent solutions of the displacement and in-plane stress components obtained using the mixed LW HSDT are in excellent agreement with those solutions obtained using the SSDQM and the BRPEM. Table 4 Table 4 that the convergent solutions are obtained when the (n x × n z ) = (256 × 8) element mesh is used.   5 show variations in the through-thickness distributions of various stress and displacement components induced in a simply supported FG beam under a sinusoidally distributed thermal load with different L/h ratios and different values of material-property gradient indices, respectively, where the LW3 theory with n l = 32 is used. In Figure 4, the length-to-thickness ratio is taken as L/h = 5, 10, and 20, and h = 0.01m, and the material-property gradient index is taken as κ e = 230. In Figure 5, they are L/h = 10, and κ e = 0, 200, and 300. A set of dimensionless variables is defined as those given in Equations (58) and (59). It can be seen in Figure 4a that the displacement component w induced in the FG beam increases when the FG beam becomes thinner, and in Figure 5a that it increases when the values of κ e become greater, which indicates the FG beam becomes softer. A comparison of the results shown in Figure 2c,d and Figure 3c,d with those in Figure 4c,d and Figure 5c,d shows that the variations in the through-thickness distributions of the transverse shear and normal stresses with the length-thickness ratio and the material-property gradient index for the thermal load cases are more drastic than the variations for the mechanical load cases. Again, the results in Figure 4c,d and Figure 5c,d show that the through-thickness distribution of the transverse shear and normal stress components τ xz and σ z appear to be a higher-order polynomial function and that the traction conditions on the top surface of the FG beam are exactly satisfied. Figures 4 and 5 show variations in the through-thickness distributions of various stress and displacement components induced in a simply supported FG beam under a sinusoidally distributed thermal load with different L/h ratios and different values of material-property gradient indices, respectively, where the LW3 theory with nl = 32 is used. In Figure 4, the length-to-thickness ratio is taken as L/h = 5, 10, and 20, and h = 0.01m, and the material-property gradient index is taken as e κ = 230. In Figure 5, they are L/h = 10, and e κ = 0, 200, and 300. A set of dimensionless variables is defined as those given in Equations (58) and (59). It can be seen in Figure 4a that the displacement component w induced in the FG beam increases when the FG beam becomes thinner, and in Figure 5a

Concluding Remarks
In this work, the authors develop the strong and weak formulations of a mixed LW HSDT for the analysis of FG beams under various boundary conditions subjected to thermo-mechanical loads, where power-law-type and exponential function-type material models are considered and the rule of mixtures is used to estimate the effective material properties of the FG beams. In the numerical examples, the strong formulation-based Navier-type analytical solutions and the weak formulationbased FEM solutions are presented. The novelty of the LW HSDT is its displacement model can capture the 3D behavior of laminated FRC and multi-layered FG beams under different boundary conditions, including the zig-zag behavior of displacement components induced in laminated FRC beams and higher-order polynomial function distributions of displacement components induced in multi-layered FG beams. The transverse shear and normal stress components induced in laminated FRC and multi-layered FG beams can be obtained using the Lagrange multipliers, which are the primary variables in the current formulation and are much more accurate than those obtained using the constitutive equations approach that has been adopted by most of the displacement-based formulations available in the literature.
Some conclusions can be summarized as follows: 1. In the LW HSDT, the transverse shear and normal stress components induced in the thermomechanically loaded FG beams can be obtained using three different approaches. In approach 1, they are obtained using constitutive equations when the displacement components are determined; in approach 2, they are obtained using stress equilibrium equations when the displacement components are determined; in approach 3, they are directly obtained using the

Concluding Remarks
In this work, the authors develop the strong and weak formulations of a mixed LW HSDT for the analysis of FG beams under various boundary conditions subjected to thermo-mechanical loads, where power-law-type and exponential function-type material models are considered and the rule of mixtures is used to estimate the effective material properties of the FG beams. In the numerical examples, the strong formulation-based Navier-type analytical solutions and the weak formulation-based FEM solutions are presented. The novelty of the LW HSDT is its displacement model can capture the 3D behavior of laminated FRC and multi-layered FG beams under different boundary conditions, including the zig-zag behavior of displacement components induced in laminated FRC beams and higher-order polynomial function distributions of displacement components induced in multi-layered FG beams. The transverse shear and normal stress components induced in laminated FRC and multi-layered FG beams can be obtained using the Lagrange multipliers, which are the primary variables in the current formulation and are much more accurate than those obtained using the constitutive equations approach that has been adopted by most of the displacement-based formulations available in the literature.
Some conclusions can be summarized as follows: 1.
In the LW HSDT, the transverse shear and normal stress components induced in the thermo-mechanically loaded FG beams can be obtained using three different approaches. In approach 1, they are obtained using constitutive equations when the displacement components are determined; in approach 2, they are obtained using stress equilibrium equations when the displacement components are determined; in approach 3, they are directly obtained using the Lagrange multipliers. Implementation of the three approaches shows that the accuracy of these approaches is in the following order: approach 2 > approach 3 > approach 1, where the symbol ">" indicates a greater degree of accuracy.

2.
The accuracy and convergence rates of the LW1, LW2, and LW3 theories are in the following order: LW3 > LW2 > LW1, where the symbol ">" indicates a faster convergence rate and a more accurate solution.

3.
Variations in the through-thickness distributions of the transverse shear and normal stresses with the length-thickness ratio and the material-property gradient index for the thermal load cases are more drastic than the variations for the mechanical load cases.

4.
Based on the convergent solutions of the LW HSDT, the accuracy of the various advanced and refined beam theories is in the following order: ESDT > SSDT > (RHSDT, TSDT) > HSDT, where the symbol ">" represents a greater degree of accuracy.