In this study, a shape optimization framework integrating reliability-based design optimization (RBDO) with isogeometric analysis (IGA) is developed. The proposed method enables shape optimization for weakly coupled thermoelastic problems by directly employing the NURBS basis functions of CAD in response analysis, while accounting for uncertainties associated with manufacturing tolerances and loads. Thermoelastic behavior frequently occurs in components operating under high-temperature and high-pressure conditions, such as reactor pressure vessels and thermal expansion pipes, where system responses are difficult to predict using empirical approaches alone. To address this issue, structural optimization is required to maximize system performance. In this work, IGA is adopted to ensure accurate geometric representation and higher continuity in the response analysis process, and solid element-based modeling is used to effectively represent complex geometries. To incorporate uncertainties arising from manufacturing and operation, reliability considerations are introduced into the design process. Since conventional reliability methods such as the first-order and second-order reliability methods (FORM and SORM) require the evaluation of design sensitivities and may exhibit reduced accuracy in highly nonlinear shape optimization problems, a sampling-based RBDO approach is employed. Surrogate models are constructed from sampled response analysis results, and reliability is evaluated using Monte Carlo simulation on the surrogate models. In the optimization framework, the control points of the IGA model are treated as probabilistic design variables to reflect manufacturing uncertainties, while temperature conditions are modeled as random variables to account for load uncertainties. The validity of the proposed approach is verified through comparison with finite element method results, and reliability-based shape optimization of a pipe example demonstrates the effectiveness and applicability of the proposed method.