The present paper considers the linear static thermal stress analysis of composite structures by means of a shell finite element with variable through-thethickness kinematic. The temperature profile along the thickness direction is calculated by solving the Fourier heat conduction equation. The refined models considered are both Equivalent Single Layer (ESL) and Layer Wise (LW) and are grouped in the Unified Formulation by Carrera (CUF). These permit the distribution of displacements, stresses along the thickness of the multilayered shell to be accurately described. The shell element has nine nodes, and the Mixed Interpolation of Tensorial Components (MITC) method is used to contrast the membrane and shear locking phenomenon. The governing equations are derived from the Principle of Virtual Displacement (PVD). Cross-ply plate, cylindrical and spherical shells with simply-supported edges and subjected to bi-sinusoidal thermal load are analyzed.Various thickness ratios and curvature ratios are considered. The results, obtained with different theories contained in the CUF, are compared with both the elasticity solutions given in the literature and the analytical solutions obtained using the CUF and the Navier’s method. Finally, plates and shells with different lamination and boundary conditions are analyzed using high-order theories in order to provide FEM benchmark solutions.