This paper describes briefly the computational algorithm that includes procedures to obtain finite-difference form of governing convection-diffusion equations, to generate an orthogonal mesh system for complex regions and to solve the finite-difference equation system, and several results of numerical simulation in comparison with experiment.The Reynolds and energy conservation equations for steady-state fully developed turbulent incompressible flows are discretized by the Efficient Finite Difference (EFD) scheme. Here secondary flow components are neglected. In the averaged energy conservation equation, anisotropic turbulent conductivity coefficients are employed based on the axial velocity distribution. An orthogonal mesh generation system has been developed that allows us to model the rod bundle geometry by assembling elementary mesh components generated for every typical sub-domain inside the flow channels. This procedure has been made efficient with a help of object-oriented programming techniques. By solving the derived equations on the boundary-fitted coordinates, good comparisons between calculation and measurement are presented in general for detailed distributions of the local shear stress, axial velocity and wall temperature in a hexagonal rod bundle in the presence of a dislocated rod. Discussion is also made on a discrepancy of the calculated wall shear stress from the experimental data near the narrowest gap position in this "geometrically disturbed" region.