Due to the difficulty in defining accurate models to represent the physical phenomena in glass furnaces. a large number of studies consider the process as a black-box system. The systemic approach proposed in this paper relies on a decomposition of the furnace into two entities, i.e. the combustion chamber and the glass bath. Next, a reactor network structure is automatically designed to reproduce the hydrodynamic features in the considered zone. These computations are based on the minimization of the quadratic difference between the network output and Residence Time Distributions (RTDs) obtained through CFD simulations. The optimization results provide a network structure showing good agreement with the RTDs computed by CFD. In addition, they highlight, in both furnace zones. some useless volumes that do not participate to the global flow. These results stand as a good basis to extend the study to the computations of kinetics, energy or environmental criteria, which in turn may be optimized. (C) 2009 Elsevier B.V. All rights reserved.