Here is a positive answer for $n=2.$ The argument doesn't seem to generalize easily to higher $n.$ The idea is to write $df\in X$ as $\star df_1 + df_2\neq 0$ and make use of Hodge decomposition.
The set $X$ consists of the matrices not of the form $(\begin{smallmatrix}a&b\\b&-a\end{smallmatrix}),$ in the $(dx,dy)$ basis (for matrices with $\sigma_1=\sigma_2,$ the middle part of an SVD is scalar, so the matrix has to be a scalar multiple of an orthogonal matrix). So we want to find an approximating sequence $g^{(n)}$ with
$$\star dg^{(n)}_1 + dg^{(n)}_2\neq 0$$
on the unit ball, with the usual Hodge star operator $\star(a\;dx + b\;dy)=-b\;dx+a\;dy.$.
For the Hodge decomposition I want to replace $\mathbb D^2$ by a more convenient space. I'll use $\mathbb R^2,$ but using a non-compact manifold is not essential an the argument could use a compactification $\mathbb R^2/\Lambda\mathbb Z^2.$
We can assume $f$ extends to a function $\bar f$ in $C^\infty_c(\mathbb R^2,\mathbb R^2),$ for example using the Sobolev extension theorem plus a standard result on density of compactly-supported smooth functions in Sobolev spaces.
The combination $\zeta=\star d\bar f_1+d\bar f_2$ is a compactly-supported smooth $L^2$ vector field.
By a perturbation we can arrange that, on a neighborhood of $\mathbb D^2,$ $\zeta$ is non-zero except at isolated points Specifically...
- pick a bounded open neighborhood $U$ of $\mathbb D^2$
- pick a $\psi$ in $C^\infty_c(\mathbb R^2,\mathbb R^2)$ that is strictly positive on $U$
- define $\phi:U\times \mathbb R^2\to\mathbb R^2$ by $\phi(x,M)=(M-\zeta(x))/\psi(x)$
- and consider a regular value $N\approx (0,0)$ for the restriction $\phi|_{U\times\{(0,0)\}}$
The preimage $\phi^{-1}(\{N\})$ is the graph $\{(x,\zeta(x)+N\psi(x))\}.$
The preimage $\phi|_{U\times\{(0,0)\}}^{-1}(\{N\})$ consists of isolated points $(x,(0,0))$ such that $\zeta(x)+N\psi(x)=(0,0).$ Projecting from the graph $\{(x,\zeta(x)+N\psi(x))\}\subset U\times \mathbb R^2$ to $U$ is a diffeo, so projecting a set of isolated points gives a set of isolated points. So the points $x\in\mathbb D^2$ with $\zeta(x)+N\psi(x)=(0,0)$ are isolated.
By pushing these out of the unit ball - we can approximate $\zeta$ in $L^2$ by a sequence of smooth $L^2$ vector fields $\gamma_n$ such that $\gamma_n\neq 0$ everywhere in the unit ball.
Each $\gamma_n$ has an orthogonal Hodge decomposition which we can write as $\gamma_n=\star dg^{(n)}_1 +dg^{(n)}_2$ where $g^{(n)}_1,g^{(n)}_2$ are determined up to additive constants.
The components $\star dg^{(n)}_1$ and $dg^{(n)}_2$ are "longitudinal and transverse" fields defined by pointwise projections in Fourier space, and since $\gamma_n$ has bounded Sobolev norms $(\int(1+|\xi|^2)^k|\hat\gamma_n(\xi)|^2d\xi)^{1/2}$ (where $\hat \cdot$ is Fourier transform), so do $\star dg^{(n)}_1$ and $dg^{(n)}_2.$ So they're smooth.
The functions $g^{(n)}_i$ are unique if we add the requirement $\int_{\mathbb D^2} g^{(n)}_i=\int_{\mathbb D^2} f_i$ for $i=1,2.$ Because Hodge decomposition is an orthogonal decomposition, $\star dg^{(n)}_1+dg^{(n)}_2\to \star d\bar f_1+d\bar f_2$ in $L^2$ implies $dg^{(n)}\to d\bar f$ in $L^2.$ The Poincaré–Wirtinger inequality then gives $g^{(n)}|_{\mathbb D^2}\to f$ in $W^{1,2}(\mathbb D^2,\mathbb R^2).$
