Unsaturated porous media constitutes a prevalent scenario in underground formations, such as oil and gas reservoirs and CO2 geological storage. The two-fluid saturated porous media is a typical manifestation of unsaturated porous media. Understanding how the elastic waves propagate in such media is fundamental for geophysical detection and security monitoring. In this work, leveraging Lo's three-phase porous media theory, a numerical simulation algorithm is introduced to analyze the propagation characteristics through such media in the time-space domain. For this purpose, the constitutive equations and the equations of motion are reformulated into velocity-stress format. By merging the time-splitting scheme with the staggered-grid finite-difference algorithm, the acoustic responses within these porous media are efficiently simulated. Moreover, reference solutions for purely dilatational point and line sources are derived to corroborate the numerical simulation results. Based on these procedures, generation mechanisms of slow compressional waves, influencing factors, and wave conversion phenomena within the interface of porous media are analyzed in focus. The study identifies four distinct waves in two-fluid saturated porous media, encompassing three types of compressional wave (P1, P2, P3) and one type of shear wave (S). Slow compressional waves (P2, P3), attributed to the relative motion in solid-fluid and fluid-fluid, respectively, are significantly affected by the fluid viscosity and frequency. Although slow waves may not be directly observed in seismic frequencies, their energy distribution needs to be emphasized, especially considering the conversion of slow waves into normally propagating P1 and S waves at interfaces.