Turbulent mixing in regular and fractal grid turbulence is investigated in this work by using direct numerical simulation (DNS). Two types of turbulence-generating grids are used: a biplane square grid (regular grid) and a square fractal grid. The thickness ratios t(r) of the fractal grids are set at 5.0 and 8.5. The grid solidity is maintained at sigma = 0.36 for all the grids. The mesh Reynolds number, Re-M = U0Meff/nu, is set at 2500 for all cases, where U-0 is the cross-sectionally averaged mean velocity; M-eff, the effective mesh size; and nu, the kinematic viscosity. The grids are numerically generated using the immersed boundary method at 4M(eff) downstream of the entrance to the computational domain. The computational domain size normalized by M-eff is 64 x 8 x 8 in the streamwise, vertical and spanwise directions for the regular grid and 64 x 16 x 16 for the fractal grids. Scalar mixing layers that initially have a step profile develop downstream of the grids. The Prandtl number is set at Pr = 0.71 considering the heat transfer in air flow. Instantaneous temperature fields, instantaneous fluctuating temperature fields and fundamental turbulent statistics are presented. The results show that turbulent mixing is more strongly enhanced in fractal grid turbulence than in regular grid turbulence for the same Re-M. In fractal grid turbulence, turbulent mixing is more strongly enhanced at t(r) = 8.5 than at t(r) = 5.0.