I'm not an expert in Vehicle Routing Problems, maybe someone else will have something more relevant to propose.
I think that a good starting point is this article:
"Efficiently solving very large-scale routing problems" (Arnold et al., 2019) DOI PDF
This paper is not exactly about your problem, but about the Capacitated Vehicle Routing Problem (CVRP), which is a simplified version of your problem.
However, they solve instances with up to 30000 locations, which is not far from your target, and, in Section 4.1, they describe how to adapt the classical Clark and Wright heuristic to large scale CVRPs and their experiments show that it can find pretty good solutions quickly (6% gap on average). In addition, it should be rather easy to implement.
Now the question is, how to adapt this heuristic to your variant? I don't see an obvious way and I am not aware of something similar in the literature about the Heterogeneous Fleet Vehicle Routing Problems.
A simple approach would be to run the heuristic with only the vehicle type with the highest capacity, and after, try to change the vehicle type of each route for a cheaper one
Note that an implementation in R will likely be slower, you might expect a factor 10 compared to the computation times from the papers
The answers about cluster-first route-second approaches reminded me of another relevant approach, route-first cluster-second, as described in this article:
"Route first—Cluster second methods for vehicle routing" (Beasley, 1983) DOI PDF
The idea is to first solve a Travelling Salesman Problem with all nodes to get a giant tour, and then to partition this giant tour into routes by dynamic programming. The author discusses how to adapt it to the case of an heterogenous fleet at the end of the article.
However, for a instance with 40000 locations, this requires to have a well optimized algorithm to solve the Travelling Salesman subproblem