We discuss novel algorithm for simulation of voxels containing multiple materials in software breast phantoms. A modification to our previous simulation of breast anatomy is proposed to improve the quality of simulated projections .  Previously, we developed a novel recursive partitioning algorithm for breast anatomy simulation in which each phantom voxel was assumed to contain a single tissue type.  As a result, phantom projection images displayed notable artifacts near the borders between regions of different materials, particularly at the skin-air boundary.  We provide theoretical evidence that  improvement in image quality without reducing voxel size is achievable by the simulation of partial volume averaging.  Thus, voxels containing more than one simulated tissue type are allowed. The linear x-ray attenuation coefficient of the voxels is the sum of the linear attenuation coefficients weighted by the voxel subvolume occupied by each tissue type. To calculate sub volumes, a local planar approximation of the boundary surface is employed. In the two-material case, the subvolumes in each voxel are computed by decomposition into simple geometric shapes.  In the three-material case, by application of the Gauss-Ostrogradsky theorem, the 3D subvolume problem is converted into one of a few simpler 2D surface area problems.  An efficient encoding scheme is proposed for the type and proportion of simulated tissues in each voxel.  Statistical methodology is introduced to validate the implementation of the algorithm and evaluate the quantitative error of the approximation.  We illustrate the proposed methodology on phantom slices and simulated mammographic projections.   Our results indicate that the proposed simulation has improved image quality by reducing quantization artifacts.