How does one solve a PDE with multiple variables, such as Stokes equation?

Question

Take one of the Stokes equations and continuity equations. Note that $\textbf{u} = (u,v,w)$ and $\textbf{f} = (f_x, f_y, f_z)$:

$\eta ({\partial^2u \over{\partial x^2}}+{\partial^2u \over{\partial y^2}}+{\partial^2u \over{\partial z^2}}) - {\partial p \over{\partial x}} + f_{x} = 0$

$\eta ({\partial^2v \over{\partial x^2}}+{\partial^2v \over{\partial y^2}}+{\partial^2v \over{\partial z^2}}) - {\partial p \over{\partial y}} + f_{y} = 0$

$\eta ({\partial^2w \over{\partial x^2}}+{\partial^2w \over{\partial y^2}}+{\partial^2w \over{\partial z^2}}) - {\partial p \over{\partial z}} + f_{z} = 0$

$ {\partial u \over{\partial x}}+{\partial v \over{\partial y}}+{\partial w \over{\partial z}} = 0$

Usually with PDEs, if we were solving for a single variable, lets say one labelled $w$, then we can solve it through something like separation of variables $w = X(x)Y(y)Z(z)$ and go from there.

However, in the case of the first three equations, we have two types of differentials; those in $u$ and those in $p$. There are two variables we would need to solve. How does one approach this?

Answer 1

Here is an outline of how to start. Take the divergence of the first three equations, and use the fourth, to get $$-\Delta p+\operatorname{div}{\bf f} =0.$$ Assuming that ${\bf f}$ is given, this no longer contains ${\bf u}$, and is the well-studied Poisson equation. How to solve it depends on the region and boundary conditions, for example Fourier transform if all of $\Bbb R^3$. Next, knowing $p$, you have Poisson equations for each of $(u,v,w)$, that still have to be linked using the continuity equation.

Answer 2

NB!! This answer briefly uses the Einstein summation convention.

---

At least in the case of the Incompressible Navier-Stokes Equations (commonly reffered to as the INSE in the business), it is possible to generate another equation that "closes" the system, though this in general is not this simple. The NSE, in most general form, is $$\begin{matrix}\rho\frac{\mathrm D\boldsymbol u}{\mathrm D t}=\nabla\cdot \boldsymbol \sigma+\boldsymbol f & \text{(momentum)} \\ \frac{\mathrm D\rho}{\mathrm Dt}+\rho\nabla\cdot\boldsymbol u=0 & \text{(continuity)}\end{matrix}$$ Here $\boldsymbol u$ is the fluid velocity, $\rho$ is the density, $\boldsymbol f$ is some external force, and $\boldsymbol \sigma$ is the Cauchy stress tensor (Google it). Under the assumptions of an incompressible liquid of constant density and in the absence of an external force, these reduce to $$\partial_t\boldsymbol u+\boldsymbol u\cdot \nabla\boldsymbol u-\nu\nabla^2\boldsymbol u=-\nabla p\tag{1}$$ $$ \nabla\cdot \boldsymbol u=0\tag{2}$$ Where $\nu=\mu/\rho$ and $p=P/\rho$ are the "kinematic" viscosity and pressure. Here the conservation of mass requirement reduces rather nicely to a requirement that $\boldsymbol u$ is divergence-free. Notice what happens when we take the divergence of $(1)$: $$\nabla\cdot\bigg(\partial_t\boldsymbol u+\boldsymbol u\cdot \nabla\boldsymbol u-\nu\nabla^2\boldsymbol u\bigg)=-\nabla\cdot\nabla p$$ The divergence commutes with both the time derivative and the Laplacian, so this becomes $$\partial_t(\nabla\cdot\boldsymbol u)+\nabla\cdot(\boldsymbol u\cdot\nabla \boldsymbol u)-\nu\nabla^2(\nabla\cdot \boldsymbol u)=-\nabla^2 p$$ The first and third terms on the LHS drop out due to $(2)$, so we get $$\nabla\cdot(\boldsymbol u\cdot\nabla \boldsymbol u)=-\nabla^2 p$$ We can deal with this rather easily if we allow ourselves to use index notation briefly: $$\nabla_i(u^j\nabla_j u^i)=-\nabla^2 p \\ (\nabla_i u^j)(\nabla_j u^i)+u^j\nabla_i\nabla_j u^i=-\nabla^2 p$$ We live in flat space, so we can commute the covariant derivatives, so $$u^j\nabla_i\nabla_j u^i=u^j\nabla_j\nabla_iu^i=(\boldsymbol u\cdot \nabla)(\nabla\cdot\boldsymbol u)= (\boldsymbol u\cdot \nabla)(0)=0$$ So the above equation reduces to $$\nabla^2 p=-\operatorname{tr}(\nabla\boldsymbol u\otimes \nabla \boldsymbol u)\tag{3}$$ Where I have defined the trace of a symmetric $(p,p)$ tensor $\mathbf T$ as $$\operatorname{tr}\mathbf T=T^{i_1\dots i_p}{}_{i_1\dots i_p}$$

---

Equation $(3)$ is known as the Pressure Poisson equation and is of massive importance in the field of computational fluid dynamics. This equation "closes" the previously "open" INSE system. We could perhaps consider writing the system in full for the INSE: $$\begin{matrix}\partial_t\boldsymbol u+\boldsymbol u\cdot \nabla\boldsymbol u-\nu\nabla^2\boldsymbol u=-\nabla p & \text{(consv. of momentum)} \\ \nabla\cdot \boldsymbol u=0 & \text{(consv. of mass)} \\ \nabla^2 p=-\operatorname{tr}(\nabla\boldsymbol u\otimes \nabla\boldsymbol u) & \text{(Pressure Poisson)}\end{matrix}$$

But since the first two equations imply the third, this is typically not done.