Multiple rounds of bipartite matching problem
You can solve this by a mild extension of the standard reduction to network flow: increase the weight of the edges touching the source and the sink to $r$. Any integer flow corresponds to a union of $r$ matchings:


An integer flow corresponds to a bipartite graph of maximum degree $r$.
We can complete this bipartite graph to an $r$-regular bipartite graph by adding vertices and edges:


Add vertices so that both sides have the same number of vertices.
If the graph is not $r$-regular, then there must be two vertices $x,y$ (one on each side) with degree less than $r$. Add $r$ new vertices $x_1,\ldots,x_r$ and $y_1,\ldots,y_r$ to each side (respectively), and connect all of them except $x_1$ and $y_1$. Connect $x$ to $y_1$ and $y$ to $x_1$. This reduces the overall deficiency of the graph $\sum_v (r - \deg(v))$ by 1.

Repeated applications of Hall's theorem shows that the new graph can be decomposed into $r$ perfect matchings.
Restricting these matchings to the original vertices and edges, we obtain a union of $r$ matchings.


Schrijver mentions this problem in his Combinatorial Optimization (Vol 1):


  Corollary 21.4b. Let $G = (V,E)$ be a bipartite graph and let $k \in \mathbb{Z}_+$. Then the maximum size of the union of $k$ matchings is equal to the minimum value of $k|V \setminus X| + |E(X)|$ taken over $X \subseteq V$.


This follows by an application of the min cut max flow theorem.