Equilibrium equations for 3D stress

Consider an infinitesimal element $dx\,dy\,dz$ in equilibrium with body forces per unit volume $X, Y, Z$ (e.g. self-weight). Stresses vary from face to face, e.g. the normal stress on the face at $x+dx$ is $\sigma_x+\dfrac{\partial\sigma_x}{\partial x}dx$.

        y
        ^     sy + dsy
        |   +----------+
  txy   |   |          | sx + dsx
  <-----|   |  dx dy   |---->
   sx   +---+----------+---> x
            sy

Equilibrium of forces in the $x$ direction

Forces = stress $\times$ area. Faces normal to $x$ have area $dy\,dz$, those normal to $y$ have area $dx\,dz$ and those normal to $z$ have area $dx\,dy$:

$$\Big(\frac{\partial\sigma_x}{\partial x}dx\Big)dy\,dz+\Big(\frac{\partial\tau_{yx}}{\partial y}dy\Big)dx\,dz+\Big(\frac{\partial\tau_{zx}}{\partial z}dz\Big)dx\,dy+X\,dx\,dy\,dz=0$$

(the terms $\sigma_x$, $\tau_{yx}$, $\tau_{zx}$ on opposite faces cancel). Dividing by $dx\,dy\,dz$:

$$\frac{\partial\sigma_x}{\partial x}+\frac{\partial\tau_{xy}}{\partial y}+\frac{\partial\tau_{xz}}{\partial z}+X=0$$

using $\tau_{yx}=\tau_{xy}$ and $\tau_{zx}=\tau_{xz}$ (from moment equilibrium).

Similarly for $y$ and $z$

$$\begin{aligned} \frac{\partial\tau_{xy}}{\partial x}+\frac{\partial\sigma_y}{\partial y}+\frac{\partial\tau_{yz}}{\partial z}+Y&=0\\ \frac{\partial\tau_{xz}}{\partial x}+\frac{\partial\tau_{yz}}{\partial y}+\frac{\partial\sigma_z}{\partial z}+Z&=0 \end{aligned}$$

Moment equilibrium about the three axes gives $\tau_{xy}=\tau_{yx}$, $\tau_{yz}=\tau_{zy}$, $\tau_{zx}=\tau_{xz}$ (symmetry of the stress tensor).

In compact form: $\ \sigma_{ij,j}+b_i=0$.