Direct flame impingement (DFI) furnaces consist of large arrays of high velocity combusting jets with temperatures up to 1700 K and impinging on complex configuration surfaces of the work pieces. This results in serious convergence problems DFI modeling and computational efforts. A new method of modeling convective-diffusion transfer (GDT) and zone radiation transfer (RT) employing different calculation schemes with a multi-scale grid is presented. Relatively coarse grid calculation domain allows use of conservative and accurate zone radiation transfer method with only modest computational efforts. A fine grid calculation domain is used to predict corrective -diffusion transfer for a representative furnace section, containing a small number of jets that allows to significantly decrease the computer time. The main difficulty of coupling between convective-diffusion transfer (CDT) and radiation heat transfer numerical computations is successfully overcome using a relatively simple algorithm. The method allows one to model the physicochemical process taking place in the DPI and reveals as well as explains many features that are difficult to evaluate from experiments. The results were obtained for high velocities (up to 400 m/s) and high firing rates. Maximum (available for natural gas-air firing) total heat fluxes up to 500 kW/m2 and corrective heat fluxes of up to 300 kW/m2 were obtained with relatively 'cold' refractory wall temperatures not exceeding 1300 K. The combustion gas temperature range was 1400-1700 K. A simplified analysis for NOx emissions has been developed as post-processing and shows extremely low NO x emissions (under 15 ppm volume) in DPI systems. Good agreement between measurements and calculations has been obtained. The model developed may be regarded as an efficient tool to compute and optimize industrial furnaces designs in limited time.