Material processing using ultra-short pulse duration lasers of low fluence produces no plasma and inflicts no damage to the material being treated. The benefit provided by the brief heating time on the order of femtoseconds also enables precision control over the spread of the heat affected zone. In this paper, the interaction between femtosecond lasers and silicon wafer is investigated. A multiscale two dimensional axisymmetric model that governs the transport dynamics in silicon is presented based on the relaxation-time approximation of the Boltzmann equation. Temperature-dependent multi-phonons, free-carrier absorptions, and the recombination and impact ionization processes are considered using a set of balance equations through which the spatial and temporal evolutions of the electron and lattice temperatures along with the electron-hole carrier density are obtained. The mechanical response of the lattice, on the other hand, is described by momentum equations. To solve the model of 17 coupled time-dependent partial differential equations without having to be concerned with non-physical oscillations in the solution, an implicit finite difference scheme on a staggered grid is developed. Unlike the conventional finite difference method in which primary variables are evaluated at grid points, the staggered finite difference scheme allows velocities and first order spatial derivative terms to be calculated at locations midway between two consecutive grid points, and shear stresses to be evaluated at the center of each element. A multi-time-scale approach involving the use of varying time steps ranging from 2.5fsec to 5psec is implemented to successfully obtain time integration results up to 10nsec. The temporal evolution and two dimensional spatial distributions of electron and lattice temperatures are presented in the first few hundreds of picoseconds when thermal equilibrium is reached. It is shown that the calculated material damage threshold based on the carrier density criterion agrees well with the published experimental data, thus validating the feasibility of the model in describing the various thermomechanical responses subject to femtosecond heating of low fluence input. The displacement and velocity distributions explain how the particles initially located at the grid nodes move as time increases and how the mechanical waves are generated by the laser heating. The thermal stress waves induced in a short time scale are proved to be highly dispersive and characteristically of broadband, low amplitude, and extremely high frequency.