To represent the formation of fuzzy nanostructures produced on a tungsten surface by exposure to a helium plasma, we have developed a hybrid simulation method that combines the binary collision approximation, molecular dynamics, and kinetic Monte Carlo calculations (BCA-MD-KMC). Since the MD code has been parallelized using the domain decomposition method (DDM) for execution in a multi-CPU environment, we developed the BCA code from scratch to mesh it efficiently with the DDM. The BCA-MD-KMC hybrid simulation code achieved a helium irradiation time of 0.1 seconds or longer, in spite of functioning at the level of atomic-scale models. In consequence, we have been able to observe the formation of concave and convex structures on a tungsten surface in the simulation.