This paper is concerned with numerical simulation of two-phase flows in complex computational regions. Both nozzle flow and jet-obstacle interaction are considered. The presence of dispersed phase (solid or liquid particles) may lead to specific thermal and erosional interaction of inertial particles with the nozzle walls and the obstacle material. The latter makes the conjugated problem much more complicated. Therefore, we consider the complete flow field in the nozzle-jet-obstacle system. The present work is a continuation of the recent study by the authors [1, 2]. A unified approach to the general problem of a two-phase nozzle-jet-obstacle flow is suggested. In this approach, both the continuous and dispersed phase behavior is calculated using the fixed rectangular grids. The solution of transient conduction equation in the solid is also carried out on rectangular grids. Both dynamics and heating/cooling of particles are calculated using the discrete-element method in Lagrangian variables. The computational model includes many mechanical effects such as collisions of particles with each other, reflection of particles from the wall surface and the feedback effect of the dispersed phase on the gas flow. The distinctive feature is the direct numerical simulation of dispersed phase dynamics, where each single real particle in the flow has its computational counterpart. All governing equations for continuous fields are solved on rectangular grids using a ghost-cell immersed boundary method. This method provides discretization of the appropriate boundary conditions via a procedure of polynomial approximation. Such approach works well for both the incompressible and compressible flows. Rectangular grids allow a straightforward implementation of high order TVD and ENO schemes for the numerical simulation of gas flows. The immersed boundary method is perfectly suited for the problems within a computational domain of varying geometry, since it doesn’t require rebuilding the grid after each boundary movement. This feature was successfully used in the numerical simulation of erosive destruction of the circular cylinder in the two-phase flow , where the mass carried away from the body resulted in moving boundaries. The current work incorporates the previous methods and algorithms into the software package allowing the numerical investigation of heterogeneous flows in more complex configurations.