Let \(n\ge 2\) , \(\Omega \subset \mathbb {R}^n\) be a bounded non-tangentially accessible domain (for short, NTA domain), and \(p(\cdot ):\mathbb {R}^n\rightarrow (0,\infty )\) a variable exponent function satisfying \(0<p_-\le p_+<\infty \) , where \(p_-:=\mathrm {ess\ inf}_{x\in \mathbb {R}^n}p(x)\) and \(p_+:=\mathrm {ess\ sup}_{x\in \mathbb {R}^n}p(x)\) . Assume that \(L_D\) is a second order divergence form elliptic operator having real-valued, bounded and measurable coefficients on \(L^2(\Omega )\) with the Dirichlet boundary condition. The main aim of this article is threefold. Firstly, the authors introduce the variable \(\textrm{BMO}\) space \(\textrm{BMO}^{p(\cdot )}(\Omega )\) and the “geometrical” variable \(\textrm{BMO}\) space \(\textrm{BMO}_z^{p(\cdot )}(\Omega )\) on \(\Omega \) , and then show that \(\textrm{BMO}^{p(\cdot )}(\Omega )=\textrm{BMO}_z^{p(\cdot )}(\Omega )\) with equivalent norms. Secondly, the authors prove the boundedness of the commutators of the Riesz transform \([\nabla L_D^{-1/2}]^k_b\) , with \(k\in \mathbb {N}\) and \(b\in \textrm{BMO}(\Omega )\) , on the weighted Lebesgue spaces \(L^p_\omega (\Omega )\) when \(p\in (1,2]\) and the variable Lebesgue space \(L^{p(\cdot )}(\Omega )\) when \(1<p_-\le p_+\le 2\) . Thirdly, the authors obtain the global gradient estimates in \(H^{p(\cdot )}_z(\Omega )\) with \(\frac{n}{n+1}<p_-\le p_+\le 1\) for the inhomogeneous Dirichlet problem of second-order divergence form elliptic equations on bounded NTA domains, where the “geometrical” variable Hardy space \(H^{p(\cdot )}_z(\Omega )\) is defined by restricting any element of the variable Hardy space \(H^{p(\cdot )}(\mathbb {R}^n)\) supported in \(\overline{\Omega }\) to \(\Omega \) , and \(\overline{\Omega }\) denotes the closure of \(\Omega \) in \(\mathbb {R}^n\) .