The Origin-Destination Cost Matrix is the best workflow in Network Analyst for processing thousands of routes at once. It could easily generate a table with 250,000 rows describing the travel time between each unique pair of nodes. However, OD Cost matrix does not export the actual route geometry between each origin and destination, as stated in the Esri help page (emphasis mine):
The closest facility and OD cost matrix solvers perform very similar analyses; the main difference, however, is in the output and the computation speed. OD cost matrix generates results more quickly but cannot return the true shapes of routes or their driving directions. It is designed to quickly solve large M x N problems and, as a result, does not internally contain the information required to generate route shapes and driving directions.
Since you are interested in the location of foot traffic based on these routes (perhaps using the Line Density tool), I think you are stuck using the Closest Facility workflow. Fortunately, all Network Analyst workflows that can be initiated from the Network Analyst toolbar also have corresponding geoprocessing tools for use in ModelBuilder of Python:
Make Closest Facility Layer geoprocessing tool. There are a cumbersome number of parameters required, but this is because there are so many options to tweak for Network Analyst workflows.
Add Locations geoprocessing tool. This is just like using the "Add Locations" window from the table of contents.
Solve geoprocessing tool. This is analogous to clicking "Solve" in the Network Analyst toolbar.
As is the case with all geoprocessing tools, you can drag and drop from the ArcToolbox pane (or the search pane) into a Model Builder model:
This model doesn't have an iterator, but it would be trivial to add one for rapid iteration through your routes. You would also need to add some tools to export your lines to individual feature classes somewhere on disk.