Here is a very explicit way of realizing a product of simplices $\Delta^r$ and $\Delta^s$ as well as a triangulation of their product with maximal simplices naturally indexed by words obtained from shuffling $a_1\cdots a_r$ with $b_1\cdots b_s$. The main thing this uses is that the type B Coxeter arrangement is a simplicial hyperplane arrangement.

The type B Coxeter arrangement is the hyperplane arrangement given by hyperplanes $x_i \pm x_j = 0$ and $x_i = 0$. Saying it is a simplicial arrangement is saying that the natural subdivision it induces on a unit sphere centered at the origin is a simplicial subdivision. Thus, we will take as our $r$-simplex the region in ${\bf R}^r$ given by $0\le a_1 \le a_2\le \cdots \le a_r$ and $\sum a_i^2 \le 1$. Likewise, we take as our $s$-simplex the region in a different ${\bf R}^s$ given by $0\le b_1\le b_2\le \cdots \le b_s$ and $\sum b_i^2 \le 1$. The product of these simplices lives in ${\bf R}^{r+s}$ and is the set of points satisfying all of these inequalities. But this is an intersection of a convex body containing the origin (namely the set of points with $\sum a_i^2 \le 1$ and $\sum b_i^2 \le 1$) and an intersection of half-spaces given by hyperplanes in the type B Coxeter arrangement for ${\bf R}^{r+s}$. But the fact that this type B Coxeter arrangement is a simplicial arrangement means that this region is a union of simplices. This triangulation into simplices is the subdivision of our region induced by the remaining hyperplanes in the type B Coxeter arrangement for ${\bf R}^{r+s}$.

In other words, the simplices are exactly the signed permutations on alphabet $a_1,\dots ,a_r,b_1,\dots ,b_s$ satisfying $0\le a_1\le \cdots \le a_r$ and $0\le b_1\le \cdots \le b_s$. But these are given by precisely the desired shuffle words. I hope this helps with what you want.

**Added later:** you also ask that the boundary of the triangulation be compatible with applying the same triangulation process to the boundaries of the original simplices. I think this construction above indeed has this property, if I understand correctly what you are asking for. To see this, it is convenient first (in terms of simplifying language) to replace the balls $\sum a_i^2 \le 1$ and $\sum b_i^2 \le 1$ above by the cubes $\max |a_i| \le 1$ and $\max |b_i| \le 1$, which we may do since the type B Coxeter complex will triangulate any full-dimensional convex body having the origin in its interior.

We are using a shuffled word such as $a_1a_2b_1a_3b_2b_3b_4$ to describe the simplex $0\le a_1 \le a_2 \le b_1 \le a_3\le b_2 \le b_3 \le b_4\le 1$. Now applying the boundary to this simplex means taking an alternating sum over all ways of turning one of these weak inequalities into an equality. If this equality is between some $a_i$ and some $b_j$, then the resulting face is an interior face of the triangulation -- because we can swap $a_i$ with $b_j$ to obtain a different shuffled word indexing a different simplex also having this face in its boundary. On the other hand, when a face results from an equality $a_i = a_{i+1}$ or $b_i = b_{i+1}$ or one from setting the first letter equal to 0 or the last letter equal to 1, then this is necessarily giving a boundary face in our triangulation of the product of simplices.

For comparison, when we take the boundary of our original simplex $\Delta^r $ (resp. $\Delta^s$), we take an alternating sum over ways of setting inequalities
to equalities $a_i = a_{i+1}$ or $a_1=0$ or $a_r=1$ (resp. $b_j = b_{j+1}$ or $b_1=1$ or $b_s=1$). Applying our triangulation procedure after doing this, thus gives shuffled words possessing exactly one such equality of the type that puts you at the boundary of the original triangulation. So *if* I understand correctly what you are asking for, I think you do get your desired compatibility.

thinkI once saw something like this in Lefschetz's book on Algebraic Topology $\endgroup$