Most of the previous studies on free-play aeroelasticity were based on the linear aerodynamic theory, which cannot consider the aerodynamic nonlinear effects. This paper proposes a framework of computational fluid dynamics/computational structural dynamics (CFD/CSD) coupling approach that can deal with the three dimensional aeroelastic problems with both the free-play nonlinearity and the aerodynamic nonlinearity. The fictitious mass method is used to construct the reduced structural equations of motion and the switching point is detected using the bisection method. The adaptive time step obtained by the bisection method is returned to the CFD solver so that both the structural and the fluid equations are integrated using the same time step. An all-movable wing with free-play at the root is considered for numerical studies. Results demonstrate the CFD/CSD coupling method can predict the stable limit cycle oscillation (LCO) effectively. The initial condition study shows that the LCO behavior is subcritical and the hysteresis response can be predicted in time domain effectively by the presented method. The viscous effect is shown to increase the LCO boundary and shift the LCO amplitude to a larger velocity in transonic regime. The LCO boundary is determined from subsonic to transonic Mach numbers. The transonic dip in LCO boundary is found by the CFD/CSD coupling method, but the equivalent linearization with doublet lattice method fails to predict this phenomenon. From the results of this study, the LCO boundary is shown to be 43.5% below the flutter boundary at most.