Concentrated solutions of short blunt-ended DNA duplexes, as short as 6 base pairs, are known to order into the nematic liquid crystal phases. This self-assembly is due to the stacking interactions between duplex terminals that promotes their aggregation into polydisperse chains with a significant persistence length. Experiments show that liquid crystal phases form above a critical volume fraction depending on the duplex length. We introduce and investigate via numerical simulations, a coarse-grained model of DNA double-helical duplexes. Each duplex is represented as an hard quasi-cylinder whose bases are decorated with two identical reactive sites. The stacking interaction between terminal sites is modeled short-range square-well potential. We compare the numerical results with predictions based on a free energy functional and find satisfactory quantitative matching of the isotropic nematic phase boundary and of the system structure. Comparison of numerical and theoretical results with experimental findings confirm that the DNA duplex self-assembly can be properly modeled via equilibrium polymerization of cylindrical particles. This insight enables us to estimate the stacking energy.