The rotated Cartesian coordinate method to remove the axial singularity of cylindrical coordinates in finite-difference schemes for elastic and viscoelastic waves

When modelling the propagation of 3D non-axisymmetric elastic and viscoelastic waves in cylindrical coordinates using the finite-difference time-domain method, a mathematical singularity occurs due to the presence of 1/r terms in the elastic and viscoelastic wave equations. For many years, this issue has been impeding the accurate numerical solution near the axis. In this work, we propose a simple but effective method for the treatment of this numerical singularity problem. By rotating the Cartesian coordinate system around the z-axis in cylindrical coordinates, the numerical singularity problems in both 2D and 3D cylindrical coordinates can be removed. This algorithm has three advantages over the conventional treatment techniques: (i) the excitation source can be directly loaded at r=0, (ii) the central difference scheme with second-order accuracy is maintained, and (iii) the stability condition at the axis is consistent with the finite-difference time-domain in Cartesian coordinates. This method is verified by several 3D numerical examples. Results show that the rotating the Cartesian coordinate method is accurate and stable at the singularity axis. The improved finite-difference time-domain algorithm is also applied to sonic logging simulations in non-axisymmetric formations and sources.