Based on this, pp 245, we go through these steps to discretize a gradient statement, namely $\nabla\phi$:
1- Gauss theorem reads,
$$
\int_V\nabla \phi dV = \oint_{\partial V}\phi dS
$$
2- Integral mean value theorem, applying on one cell (i.e. $c$) results,
$$
{\bar{\nabla\phi}}_c =\frac 1 {V_c} \oint_{\partial V_c}\phi dS
$$
3- Applying $2^{nd}$ order accurate discretisation for RHS, we end up with ($f$ stands for cell face centroids),
$$
{\bar{\nabla\phi}}_c =\frac 1 {V_c}\Sigma \phi_f S_f
$$
The question is :
Does applying Integral mean value theorem, namely ${\bar{\nabla\phi}}_c =\frac 1 {V_c} \int_{V_c}\nabla \phi dV$, maintain $2^{nd}$ order accuracy of discretisation?
1 Answer
Not strictly no. They will reconstruct only linear functions on simple meshes exactly. You can get a better approximation (for non stretched meshes) using weighted least squares gradient reconstruction, but that still wont be exactly second order accurate. You can find more information on the gradient reconstruction techniques and their pitfalls here:
https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20040070704.pdf