This paper addresses smooth constrained optimization problems with unknown objective and constraint functions The main goal of this work is to generate a sequence of feasible (herein, referred to as safe) primal-dual pairs converging towards a KKT pair. Assuming to have prior knowledge on the smoothness of the unknown functions, we propose a novel zeroth-order method that iteratively computes quadratic approximations of the constraint functions, constructs local feasible sets, and optimizes over them. We prove that this method returns an eta-KKT pair within O(d/eta 2) iterations and O(d2/eta 2) samples (where d is the problem dimension) while every sample is within the feasible set. Moreover, we numerically show that our method can achieve fast convergence compared with some state-of-the-art zeroth-order safe approaches. The effectiveness of the proposed approach is also illustrated by applying it to a nonconvex optimization problem in optimal control. (c) 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).