In the near future, gravitational wave detection is set to become an important observational tool for astrophysics. It will provide us with an excellent means to distinguish different gravitational theories. In effective form, many gravitational theories can be cast into an f (R) theory. In this article, we study the dynamics and gravitational waveform of an equal-mass binary black hole system in f (R) theory. We reduce the equations of motion in f (R) theory to the Einstein-Klein-Gordon coupled equations. In this form, it is straightforward to modify our existing numerical relativistic codes to simulate binary black hole mergers in f (R) theory. We considered binary black holes surrounded by a shell of scalar field. We solve the initial data numerically using the Olliptic code. The evolution part is calculated using the extended AMSS-NCKU code. Both codes were updated and tested to solve the problem of binary black holes in f (R) theory. Our results show that the binary black hole dynamics in f (R) theory is more complex than in general relativity. In particular, the trajectory and merger time are strongly affected. Via the gravitational wave, it is possible to constrain the quadratic part parameter of f (R) theory in the range |a2| < 10 11 m 2 . In principle, a gravitational wave detector can distinguish between a merger of binary black hole in f (R) theory and the respective merger in general relativity. Moreover, it is possible to use gravitational wave detection to distinguish between f (R) theory and a non self-interacting scalar field model in general relativity.