Why is the pressure gradient zero at a wall?

This usually only applies to a wall bounded flow and is normally restricted to incompressible fluids. This result usually manifests in boundary layer theory and can be obtained through order of magnitude analysis of the Navier-Stokes equations. The steady, incompressible, and constant property momentum equation in the $y$ direction takes the form, $$ u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right)$$ The order of magnitude of each term inside the boundary layer goes as follows, $$ O\left[u \frac{\partial v}{\partial x}\right] = O\left[\frac{\delta}{L^2} U_e^2\right] $$ $$ O\left[v \frac{\partial v}{\partial y}\right] = O\left[\frac{\delta}{L^2} U_e^2\right] $$ $$ O\left[\frac{1}{\rho} \frac{\partial p}{\partial y}\right] = O\left[\frac{\delta}{L^2} U_e^2\right] \text{(at most)}$$ $$ O\left[\nu \frac{\partial^2 v}{\partial x^2}\right] = O\left[\frac{\delta^2}{L}\frac{\delta}{L^2} U_e^2\right] $$ $$ O\left[\nu \frac{\partial^2 v}{\partial y^2}\right] = O\left[\frac{\delta}{L^2} U_e^2\right] $$

Where $\delta$ is the boundary layer height, $L$ is the characteristic length of the body, and $U_e$ is the external flow velocity at the edge of the boundary layer. An additional constraint is that $\delta/L \ll 1$. Notice that each term has an order of magnitude comparable to $(\delta/L^2) U_e^2$, except for the normal pressure gradient term and the viscous term from the $x$ direction. We first recognize that $\delta^2/L$ is a very small quantity and essentially removes the $\nu \partial^2 v/ \partial x^2$ term from the equation. Similarly, only at the edge of the boundary layer where the viscous forces become negligible (i.e. high Reynolds number) does the pressure gradient term order of magnitude approach $(\delta/L^2) U_e^2$. This was first observed by Prandtl, to which he deduced across the boundary layer we can write, $$ \frac{1}{\rho} \frac{\partial p}{\partial y} \approx 0$$ or more conventionally, $$ \frac{\partial p}{\partial y} \approx 0 $$

As for your second question, this only applies to the static pressure. Also, all of this assumes the flow is attached to the wall.

It comes from the notion of the boundary layer and whether it stays attached to the wall or not. If you consider the momentum equation normal to the wall, the only way there can be a pressure gradient normal to the wall is if there is a velocity or acceleration normal to the wall. If that is the case, then the boundary layer is no longer attached.

This holds as you get infinitesimally close to the wall (well, so long as it is still a continuum). At some point, the flow is attached to the wall, even if it is a tiny thin layer, and so in simulations, the gradient is zero at the wall. If you have a proper grid resolution, all is well.

Given some of the comments to this question I thought it might be useful to provide a somewhat more detailed answer. I start with a basic review of some of the mathematics, and will close with comments on approaches to the numerical solution of the problem.

Remarks on the Mathematics of the incompressible Navier-Stokes problem

To simplify matters I will consider a scenario where we restrict ourselves to two-dimensional incompressible Navier-Stokes flow, described in a Cartesian $x_1$-$x_2$ coordinate system. The velocity and pressure fields are $\mathbf u={\mathbf u}(\mathbf x,t)=(u_1,u_2)^T(x_1,x_2,t)$ and $p(\mathbf x,t)=p(x_1,x_2,t)$.

I will also assume that our problem has Dirichlet boundary conditions for the velocities everywhere in a square domain $\Omega=\{(x_1,x_2)\,\,\vert \,\,0\le x_1\le1,\, 0\le x_2\le1\}$ with boundary $\Gamma$.

The system of equations describing incompressible flow in this domain is given by

$${\mathbf\nabla}\cdot{\mathbf u}=0,$$ $$\frac{\partial\mathbf u}{\partial t} +{\mathbf u}\cdot{\mathbf\nabla}{\mathbf u}=-{\mathbf\nabla}p+\frac{1}{Re}\Delta{\mathbf u}$$

This system needs to be complemented by initial and boundary conditions, for which we choose

$${\mathbf u}(x_1,x_2,0)={\mathbf u}_0(x_1,x_2),\quad (x_1,x_2)\in\Omega,$$

and

$${\mathbf u}(x_1,x_2,t)={\mathbf u}_\Gamma(x_1,x_2,t),\quad (x_1,x_2)\in\Gamma,$$

where $\Gamma=\partial\Omega$ is the domain boundary. We allow for arbitrary ${\mathbf u}_\Gamma(x,y,t)$, with the exception of the bottom piece of the boundary, which we choose as a solid no-slip wall, so we have

$${\mathbf u}_\Gamma(x_1,x_2,t)={\mathbf 0},\quad (x_1,x_2)\in\{(x,0)\,\,\vert \,\,0\le x\le1\}.$$

Notice, first, that the above set of equations does not include, and does not require any boundary conditions for pressure. Indeed, pressure appears as a Lagrange multiplier in the momentum equations which is used to project the solution of the momentum equations onto the space of divergence-free vector fields. In other words, its function is to ensure that the continuity equation is satisfied, and mass is conserved.

It is instructive to take a closer look at the properties of the resulting pressure field. Taking the divergence of the momentum equations and rearranging, we find that the pressure field satisfies the Poisson equation

$$\Delta p=\nabla\mathbf u:\nabla\mathbf u,$$

where the colon stands for the inner product $\frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i}$ where the summation convention was used.

To find the boundary values that are assumed by this pressure field, we re-write the momentum equations at the boundaries. Since we are specifically interested in the wall-normal pressure gradient on a solid wall, let us project the momentum equations on the wall-normal coordinate $x_2$ to obtain

$$\frac{\partial p}{\partial x_2}=-\frac{\partial u_2}{\partial t} -{\mathbf u}\cdot{\mathbf\nabla}{u_2}+\frac{1}{Re}\Delta u_2.$$

Using the boundary condition $\mathbf u\equiv0$ at the wall this simplifies to

$$\frac{\partial p}{\partial x_2}=\frac{1}{Re}\frac{\partial^2 u_2}{\partial x_2^2}.$$

Thus, it is clear that the wall-normal pressure gradient is non-zero in general, and prescribing zero pressure gradient at the wall is mathematically inconsistent with the original problem. Since this point came up in a comment, I will also note that we cannot, in general, take comfort in the appearance of the $1/Re$ factor in the above expression, trying to argue that the right-hand side can be neglected for large Reynolds number. This is not so, since we are considering the near-wall region of the flow, and in order to assess the relative magnitude of said right-hand side, we may need to consider a local Reynolds number, appropriate to the size of important near-wall flow structures. Long story short, it turns out that for a vast class of important flows, such as near-wall turbulence, this local Reynolds number is of the order of unity. In other words, we cannot neglect the right-hand side for any such flow.

Let us now shine a little more light on the potential consequences of providing incorrect boundary conditions for the pressure field. This is most easily seen by considering an alternative but equivalent formulation of the Navier-Stokes problem above. It will turn out that this formulation, known as the "pressure-Poisson equation formulation" is also important for numerical approaches to the solution of these equations. For this formulation we replace the combination of continuity and momentum equations in the so-called "primitive variables formulation" of the Navier-Stokes equations by a combination of pressure-Poisson and momentum equations, like so:

$$\Delta p=\nabla\mathbf u:\nabla\mathbf u,$$ $$\frac{\partial\mathbf u}{\partial t} +{\mathbf u}\cdot{\mathbf\nabla}{\mathbf u}=-{\mathbf\nabla}p+\frac{1}{Re}\Delta{\mathbf u}.$$

This system of partial differential equations now requires the following initial and boundary conditions:

$${\mathbf u}(x_1,x_2,0)={\mathbf u}_0(x_1,x_2),\quad (x_1,x_2)\in\Omega,$$ $$\mathbf\nabla\cdot{\mathbf u}(x_1,x_2,t)=0,\quad (x_1,x_2)\in\Gamma,$$ $${\mathbf u}(x_1,x_2,t)={\mathbf u}_\Gamma(x_1,x_2,t),\quad (x_1,x_2)\in\Gamma,$$

Notice that we now need to require an additional boundary condition, since the pressure-Poisson equation formulation was obtained from the original system via a differentiation operation, thus increasing the order of the system. You can also see that the system of PDEs we have now obtained has a non-standard structure, in that we now have three boundary conditions for the velocity field, and none for the pressure. I will make two important remarks:

1. This situation leads to significant difficulties in the implementation of proper numerical schemes for the solution of this set of equations. Such schemes are possible, however, and are known as "influence matrix methods".
2. There is a shockingly large number of references in the literature who recommend using a pressure gradient boundary condition to replace the first of the boundary conditions we have stated above. There are two variants of these: (a) Use a zero pressure gradient condition. This results in a well-posed mathematical problem, but the solution of this problem cannot satisfy the continuity equation for incompressible flow. Simply put, the resulting solution is wrong. We emphasize that the solution is mathematically wrong, and a numerical scheme using this approach will not converge to the correct Navier-Stokes solution no matter what. (b) Use $\nabla p=(1/Re) \Delta\mathbf u$ as a boundary condition. This makes matters worse: The resulting problem is now mathematically ill-posed, and does not have a definite solution at all.

Finally, to better understand the role that the pressure plays in incompressible Navier-Stokes solutions, it may be instructive to look at the behavior of the divergence of the velocity field, given the pressure Poisson equation. Taking the divergence of the momentum equations and substituting the pressure Poisson equation in there, we arrive at

$$\frac{\partial}{\partial t}(\mathbf\nabla\cdot\mathbf u)=\Delta(\mathbf\nabla\cdot\mathbf u).$$

Thus the velocity divergence satisfies a heat equation. If proper boundary conditions are not enforced, based on the appropriate maximum/minimum theorems for such equations we therefore expect maximum divergence errors at the boundaries, which will diffuse away from the boundaries into the domain. It should go without saying that for many boundary layer flows (transitional and turbulent boundary layers in particular), such errors can be disastrous and completely destroy the representation of the physics of the flow.

Remarks on numerical methods

When considering approaches to the numerical solution of the Navier-Stokes problem, it is crucial to understand the exact details of the methods considered. When doing so one will find that the concept of "boundary condition" is often handled in an exceedingly sloppy manner in the community, which can make it difficult to tell what equations are in fact being used, and why. This is particularly true in the case of finite-difference and finite-volume codes. Examples of schemes that are either over-determined or mathematically ill-posed abound in this area. Spectral methods and finite-element methods typically enforce somewhat more mathematical discipline, and therefore seem to be a bit less prone to some of the pitfalls. However, there is one statement we can make with certainty:

There is no consistent and accurate numerical scheme that violates the mathematics outlined above. Incorrect boundary conditions for the pressure result in incorrect solutions, end of story.

Specifically, no successful numerical scheme in the above sense uses a zero pressure gradient boundary condition for the pressure. However, there are a number of numerical methods that look like they are using improper boundary conditions, and among these are some of the most popular schemes used in both research and commercial codes. For example, in fractional step methods, as well as the SIMPLE algorithm, iterative procedures are employed that use a zero gradient boundary condition for a pressure correction step. I want to emphasize that such methods can be perfectly legitimate, consistent with the exact equations, and converge to accurate solutions. There's a large number of possible variants of such approaches and it does not make sense to discuss these in any detail here. Suffice it to say that quite often one finds that the "pressure" used in the core of such algorithms does not represent the mathematically correct pressure field but instead is a pseudo-pressure. The resulting solution produces a consistent approximation of the velocity field, but the pseudo-pressure systematically deviates from the correct pressure field. This situation may or may not be corrected in additional steps of the algorithm.

This situation is further complicated by the fact that the precise implementation details of many numerical methods remain unclear, with the result that algorithms that simply should not work based on their formal as-published description do give useful approximations, for reasons that remain unclear, however. I will say that I am intimately familiar with the development process of such codes from my own research group and from colleagues. All too often what really happens is that the graduate student starts with some ill-conceived code that should not work, and then keeps working on it until "the code is fixed" and produces acceptable results in some sense. A discussion of what is wrong with this picture truly would lead me too far afield here...

Finally, for a much more detailed review of this topic see the review article by Rempfer in Applied Mechanics Reviews ("On Boundary Conditions for Incompressible Navier-Stokes Problems", AMR (59) 2006).