I think the closest to an "analytical" answer (and I do not mean in terms of the accuracy, but in the sense that you can solve it by using pen and paper instead of a computer) would be to use linearization. Consider $g(X_1,X_2)=e^{X_1}+e^{X_2}$, we'd now like to compute $VaR_\alpha(g(X_1,X_2))$.
Linearization of $g$ gives:$$g(X_1,X_2)\approx g(\mathbf{\mu}) + \nabla g^T(\mu)(\mathbf{X}-\mathbf{\mu})$$
Thus we have that $$VaR_\alpha(g(X_1,X_2)) \approx VaR_\alpha(g(\mathbf{\mu}) + \nabla g^T(\mu)(\mathbf{X}-\mathbf{\mu})) = VaR_\alpha(\nabla g^T(\mu)(\mathbf{X}-\mathbf{\mu})) - g(\mathbf{\mu})$$
Where the last equality holds due to the translation invariance of value-at-risk. Now we have that $$\nabla g^T(\mu)(\mathbf{X}-\mathbf{\mu})= e^{\mu_1}(X_1-\mu_1) + e^{\mu_2}(X_2-\mu_2) = e^{\mu_1}X_1+e^{\mu_2}X_2 -\mu_1e^{\mu_1} - \mu_2e^{\mu_2}$$
Once again we can move out the constants due to the translation invariance of value-at-risk. (Note that I skip the risk-free rate).
Now we have $$VaR_\alpha(g(X_1,X_2)) \approx VaR_\alpha(e^{\mu_1}X_1+e^{\mu_2}X_2)+e^{\mu_1} + e^{\mu_2} -\mu_1e^{\mu_1} - \mu_2e^{\mu_2}$$
Since $X_1,X_2$ are just normally distributed, so will their sum be. Thereafter its just straightforward computations.
Worth noting is that linear approximation only works well when the probability mass in centered closely around the mean, i.e. if tails are light and variance is small. It becomes particularly worrisome since in the value at risk we're measuring probabilities far out in the tail.