Abstract: In this paper, we proposed a novel numerical method, zonal free element method (ZFLM), and used the proposed method to compute thermal stress in composite structures. ZFLM is a strong-form numerical method which solves the governing equations in differential form. For each node, we use two (two-dimension problems) or three (three-dimensions problems) lines to form a cross-line system. Then, we use Lagrange interpolating method to interpolate nodal coordinates and approximate the variables on each line. The gradients in the curvature direction are computed by the gradients of interpolating functions along the line. By a recursive procedure, the second or higher order of derivatives can be obtained by the expressions of the first order derivatives. Substituting the expressions of derivatives into the governing partial difference equations, we obtain the discretized linear system of equations. To solve the problem involving multiple kinds of composite structures efficiently, we use a zonal method. In the zonal method, we divide the computational domain into several regular zones by material types and geometric characteristics. We insert nodes in each zone by interpolating functions and use the finite line method to assemble the discretized governing equations at these nodes. For the nodes at the interfaces which are shared by two or more zones, the traction-equilibrium equations and the compatibility conditions of variables are used to construct the linear algebraic equations. For the irregular geometries and the nodes where the loads jump, we add up the traction-equilibrium equations of each neighbor faces in different directions to improve the robustness of the proposed method. We use the proposed method to solve several thermal stress problems in two- and three-dimension. The results of test cases indicate that the proposed method has good accuracy and a significant priority in problems involving stress concentration. Because the collocation method is used, the stress on the boundary is more accurate.