In my own search for a rigourous proof I found that it is much easier to use an equivalent condition for what constitutes a simplicial complex; namely
A set of simplicies $\mathcal{S}$ in some euclidean space is a simplicial complex if and only if
- The faces of simplicies in $\mathcal{S}$ are also in $\mathcal{S}$, and
- The interiors of distinct elements of $\mathcal{S}$ are disjoint.
The fact that the face of an element of the barycentric subdivision is also an element of the barycentric subdivision is practically trivial so I will roughly lay out a proof of condition 2.
So, let $\text{Bary}(\mathcal{S})$ be the barycentric subdivision and fix $\upsilon,\rho \in \text{Bary}(\mathcal{S})$ such that $\text{int}(\upsilon) \cap \text{int}(\rho) \neq \emptyset$. We wish to show that $\upsilon = \rho$ which proves condition 2 and that $\text{Bary}(\mathcal{S})$ is a simplicial complex.
Take $\sigma_0 \subset \dots \subset \sigma_m \in \mathcal{S}$ and $\tau_0 \subset \dots \subset \tau_n \in \mathcal{S}$ such that $\upsilon = \text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_m}\})$ and $\rho = \text{Convex}(\{b_{\tau_0},\dots,b_{\tau_n}\})$.
Proceeding by induction assume that for some $k \in \{0,\dots,\min\{m,n\}-1\}$ we have $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-k}}\})) \cap \text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-k}}\})) \neq \emptyset$ and $\sigma_{m-i} = \tau_{n-i}$ for all $i \in \{0,\dots,k-1\}$.
The case $k=0$ is trivial. It is not too difficult to see, by explicitly writing out barycentric coordinates and rearranging, that $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-k}}\})) \subseteq \text{int}(\sigma_{m-k})$ and $\text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-k}}\})) \subseteq \text{int}(\tau_{n-k})$. So, by assumption and since $\mathcal{S}$ is a simplicial complex, we have $\sigma_{m-k} = \tau_{n-k}$.
Next, fix $p \in \text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-k}}\})) \cap \text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-k}}\}))$. If we take the ray from $b_{\sigma_{m-k}} = b_{\tau_{n-k}}$ through $p$ then on the one hand we find it intersects $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-k-1}}\}))$ and on the other it intersects $\text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-k-1}}\}))$. Again, one can write out $p$ in barycentric coordinates to explicitly find the intersection points.
However, we have $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-k-1}}\})) \subseteq \sigma_{m-k-1}$ which is a proper face of $\sigma_{m-k}$ and likewise $\text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-k-1}}\})) \subseteq \tau_{n-k-1}$ which is a proper face of $\tau_{n-k}$. But, $\sigma_{m-k} = \tau_{n-k}$ and, since this is a convex set, the ray from the interior point $b_{\sigma_{m-k}} = b_{\tau_{n-k}}$ through $p$ must intersect the boundary in at most one point.
Hence, the intersections in $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-(k+1)}}\}))$ and $\text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-(k+1)}}\}))$ are in fact the same so that $\text{int}(\text{Convex}(\{b_{\sigma_0},\dots,b_{\sigma_{m-(k+1)}}\})) \cap \text{int}(\text{Convex}(\{b_{\tau_0},\dots,b_{\tau_{n-(k+1)}}\})) \neq \emptyset$.
As a result we have by induction that $\sigma_{m-i} = \tau_{n-i}$ for all $i \in \{0,\dots,\min\{m,n\}\}$. If $m < n$ then we would have $\{b_{\sigma_0},\dots,b_{\sigma_m}\} \subset \{b_{\tau_0},\dots,b_{\tau_n}\}$ so that $\upsilon$ is a proper face of $\rho$ which contradicts the fact that $\text{int}(\upsilon) \cap \text{int}(\rho) \neq \emptyset$. Likewise, we cannot have $n < m$.
As such, $m = n$ so that $\sigma_i = \tau_i$ and $b_{\sigma_i} = b_{\tau_i}$ for all $i \in \{0,\dots,n\}$. Thus, $\upsilon = \rho$ and therefore $\text{Bary}(\mathcal{S})$ is a simplicial complex.