Computational micromechanics, aimed at performing numerical simulations at the micro scale, is of utmost importance, particularly in the current context of digital twins in the mechanics of materials. The objective in this work is to develop an efficient and accurate numerical method capable of simulating geometrically complex composite microstructures with numerous inclusions connected to the matrix through various non-linear interface behaviors. To achieve this goal, we employ IsoGeometric Analysis (IGA), which provides increased per-degree-of-freedom accuracy, and leverage the recently introduced immersed boundary-conformal method. This method allows us to bypass the challenging meshing task while ensuring conformal matrix-to-inclusion interfaces through the construction of conformal layers from the latter. Hence, we separate the difficulties by creating (i) non-linear but conformal interfaces between the matrix and the inclusions, and (ii) non-conformal but perfect interfaces to couple the introduced layers within the matrix and inclusions. Then, the approach is enhanced with the Large Time Incremental (LaTIn) method, enabling the separation of non-linear interface equations from those related to the subdomains. The latter are linear and subdomain-wise independent, facilitating parallelized solutions. LaTIn solvers have previously been applied, e.g., with FE immersed approach for multiple contacts. This work extends LaTIn to higher-order immersed IGA. The non-conformal coupling is eventually addressed using Nitsche’s approach, which has demonstrated optimal convergence behavior for splines. The accuracy and efficiency of the developed algorithm are demonstrated by solving a range of non-linear examples in 2D, including varying numbers of inclusions in unilateral contact, frictional contact, and delamination within the matrix of the composite microstructure.