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)\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 concrete, let $X^c$ denote $\{(\begin{smallmatrix}a&b\\b&-a\end{smallmatrix})\mid a,b\in\mathbb R\},$ choose a compactly supported smooth function $\psi$ that is positive on a bounded open neighborhood $U$ of $\mathbb D^2,$ and pick a small regular point $N$ for the smooth map $\phi:U\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.$ Projecting isolated points from the graph $\{(x,\zeta(x)+N\psi(x))\}$ to $U$ gives a set of isolated points. So the points $x\in\mathbb D^2$ 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 $\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).$