Metastable γ′′ (D022-Ni3Nb) particles, exhibiting three variations in orientation, serve as strengthening precipitates in the commercial Inconel 625 Ni-based superalloy. Understanding the morphological evolution and fractions of these particles within a Representative Volume Element (RVE) matrix is crucial for evaluating the hardening effects and guiding the design of additive manufacturing processes to improve yield strength. A computational coupled phase-field model was developed to predict solid-state phase transformation kinetics within mechanical parts during the heat-treating process. The time evolution of the concentration fields during precipitation is governed by the generalized diffusion equation, i.e., the Cahn–Hilliard equations. Simultaneously, the time evolution of the structural order parameter (phases) is governed by the time-dependent Allen-Cahn equation. The precipitation model is quantitative, utilizing ab initio calculations of elastic constants, experimental data on lattice parameters, precipitate matrix orientation relationships, interfacial energy of each individual precipitate phase, and inter-diffusivities. Additionally, a Ni–Nb–Al pseudo-ternary thermodynamic database specifically developed for IN625 is employed. Simulation results illustrate how alloy composition, lattice misfit, external stress, temperature, and time influence precipitate microstructure and variant selection during isothermal heat history fluctuations. Importantly, the model makes no a priori assumptions about key microstructural features, including size, shape, volume fraction, and spatial distribution of different types of precipitates and different variants of the same precipitate phase.