Chemotaxis-driven aggregation and transport in biological media are frequently accompanied by fluid motion, and image-derived observations motivate reliable coupled PDE solvers for quantitative biomechanics modeling. We consider a thermodynamically consistent Keller--Segel--Navier--Stokes (KSNS) system derived within the Energetic Variational Approach (EnVarA), where the chemotaxis subsystem involves a logarithmic free-energy potential that is nonlinear and singular as the cell density approaches vacuum. Standard discretizations may violate positivity and generate nonphysical energy growth, leading to severe instability in long-time simulations. We propose a second-order accurate, structure-preserving numerical scheme for the KSNS system. The Keller--Segel part is reformulated as a coupled $H^{-1}$ gradient flow with non constant mobility and an $L^{2}$ gradient flow, and is discretized in time by a modified Crank--Nicolson method. Carefully designed artificial regularization terms enforce positivity of the cell density at the fully discrete level. The incompressible Navier--Stokes equations are advanced by a second-order semi-implicit temporal discretization and a marker-and-cell (MAC) staggered-grid spatial discretization, yielding a discretely divergence-free velocity field and stable chemo-fluid coupling. Rigorous analysis proves unique solvability and preserves three key structural properties: positivity of cell density, exact mass conservation, and monotone dissipation of the total free energy. Using higher-order asymptotic expansions with rough and refined error estimates, we establish optimal convergence rates. Numerical experiments confirm robust long-time behavior, correct energy decay, and optimal accuracy under mesh-time refinement.