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.$.
The Hodge decomposition I want to replace $\mathbb D^n$ by a more convenient space. I'll use $\mathbb R^n,$ but using a non-compact manifold is not essential an the argument could use a compactification $\mathbb R^n/\Lambda\mathbb Z^n.$
We can assume $f$ extends to a function $\bar f$ in $C^\infty_c(\mathbb R^2,\mathbb R^2)\cap W^{1,2}(\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 assume that $\zeta$ is non-zero except at isolated points. To be specific, let $X^c$ denote $(\begin{smallmatrix}a&b\\b&-a\end{smallmatrix}),$ choose a compactly supported function $\psi$ that is positive on $\mathbb D^2,$ and pick a sequence of regular points $N$ for the smooth map $\phi:\mathbb R^n\times X^c\to \mathbb R^{2\times 2}$ defined by $\phi(x,M)=(M-\zeta(x))/\psi(x).$ The preimage consists of isolated points $(x,M)$ such that $\zeta(x)+N\psi(x)=M\in X^c.$ The set of $x$ values in this preimage that fall in $\mathbb D^n$ is compact - the possible $M$ are bounded because $\zeta$ is bounded. So the points $x\in\mathbb D^n$ with $\zeta(x)+N\psi(x)\in\ X^c$ are isolated.
By pushing these out of the unit ball - composing with a suitable smooth diffeo $\mathbb R^2\to\mathbb R^2$ that affects only a region of small measure - 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 the Fourier transform of $\gamma_n$ decays super-polynomially, so does the Fourier transform of $\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).$