We consider the approximation of two coupled wave equations with internal damping. \[ \begin{cases} u_{tt}-u_{xx}+\alpha y+\beta u_t=0&\text{in }(0,1)\times\mathbb R_+,\\ y_{tt}-y_{xx}+\alpha u+\gamma y_t=0&\text{in }(0,1)\times\mathbb R_+,\\ u(0,t)=u(1,t)=y(0,t)=y(1,t)=0,&\text{for all }t>0,\\ u(\cdot,0)=u_0,\;u_t(\cdot,0)=u_1,\;y(\cdot,0)=y_0,\;y_t(\cdot,0)=y_1&\text{in }(0,1). \end{cases} \] Our goal is to damp the spurious high frequency modes by introducing numerical viscosity terms in the approximation scheme. With these viscosity terms, we show the exponential or polynomial decay of the discrete scheme when the continuous problem has such a decay (since the spectrum of the spatial operator associated with the undamped system satisfies the generalized gap condition).