Let $X \sim \mathcal{N}(\mu,\sigma^2)$. I think it's true that

$$\mathbb E \Phi(X) = \Phi \left(\frac{\mu}{\sqrt{1+\sigma^2}}\right)$$ where $\Phi$ is the cdf of standard normal.

This holds up under simulation: if you run the R code

mu <- rnorm(1)
sigma <- abs(rnorm(1))
number_of_simulations <- 1e5
simulated_result <- mean(pnorm(rnorm(number_of_simulations,mu,sigma)))
analytical_result <- pnorm(mu*sqrt(1/(1+sigma^2)))
abs(simulated_result - analytical_result)

you should get a number close to 0.

Let $f$ be the pdf of $X$. We have $$\mathbb E \Phi(X) = \int_{-\infty}^{\infty}\Phi(x)f(x)\,dx$$ which doesn't look trivial to solve. Is there a proof or a well known result that can prove this? Or is there a good intuitive explanation of it?

The question originally arose from the Merton framework for credit risk. Define $$Q = \frac{B - \sqrt{\rho}Z}{\sqrt{1-\rho}}$$ where $Z\sim \mathcal{N}(0,1)$ then we get $\mathbb E\Phi (Q) = \Phi(B)$ which is equivalent to the formula in this question.