Heterogeneity in engineering materials arises from various factors spanning material science and manufacturing [1], such as reinforcing inclusions, defects like inhomogeneities and voids, necessitating rigorous investigation [2]. Finite Element Analysis (FEA) captures multiscale heterogeneous material behaviors but faces accuracy-efficiency tradeoffs and convergence challenges from the stochastic distribution of inhomogeneities and contact nonlinearity, with solution fidelity reliant on mesh discretization quality and computational resources. Experimental techniques provide microstructural insights but lack resolution for heterogeneous responses under contact loading. This work presents a robust and efficient approach to addressing material nonlinearity using an inclusion-based framework. The Numerical Equivalent Inclusion Method (NEIM) is employed to model arbitrary inhomogeneities in non-homogeneous half-space. A single-loop iteration schemed based on the conjugate gradient method (CGM) is utilized to solve the contact problem with various inhomogeneities. The numerical integration of NEIM with the Fast Fourier Transform (FFT) algorithms leverages computational advancements, ensuring rapid and reliable solutions of contact behaviors in layered heterogeneous materials. Fundamentally, the combination of NEIM and FFT principles in contact problems not only enhances computational accuracy and efficiency but also offers a versatile and powerful solution. This approach represents a significant advancement in tackling the complexities of these problems, providing valuable insights for engineers and researchers across various disciplines.