You most probably do not need an answer anymore. However I stumbled in just the same kind of inconsistency and I found this post, so this comment might help others like me. I understand this is the configuration you have in mind... very commonly used to introduce shear stress:
You are right: on the one hand, you can't have any force on a free surface (like the vertical ones on the left and right side of the cube) in a stationary configuration; on the other one you need a symmetric stress tensor otherwise you get problems with the torque of surface forces... no doubts about that... but in this way you also get a non-zero force on those vertical surfaces. Good point, it does not make any sense: I agree.
Well, today I struggled for hours trying to understand what I was not getting here and... Boy - I was kind of shocked in the beginning - the problem is in a stupid assumption, i.e. that what is show in the figure above is REALLY what happens if you pull a cube laterally from the top surface and at the same time prevent the bottom surface from sliding on the floor. The problem is that you simply won't get that nice uniform diamond-like deformation... but something quite more complicate and definitely non-uniform. So in the end the reason of the inconsistency you mention is that the "solution" in the figure is not a solution of the mechanical problem.
In particular (beam theory can help here) I would say that:
If you consider the bottom surface, it really has to be glued to the ground because the left edge will be pulled up and the right one will be pushed against the floor. Lateral friction is not enough to prevent the cube from rotating: try with a foam cube and you will see...
The bottom surface will exert a complex force profile with various components perpendicular to the ground (both positive and negative depending on the position). THIS will balance the torque of the force on the top surface if external forces just pull the top surface of the cube.
Obviously the strain is not at all uniform here... there must be some sort of non-trivial bending, in particular at the base of the cube.
Lateral surfaces are free and, you are right, they have to exert zero force! Be sure that force on those faces is indeed zero in the correct solution of this problem.
In short, once you properly glue the bottom of the cube to the floor, you can expect something more similar to what visible here
rather than in the picture above.
It is just funny that sketches like the one above are so widely used to explain shear stress. I think it is kind of risky and can lead to inconsistencies and (legitimate) doubts in students, I wonder how many persons are aware.
I have to say that some texts indicate the presence of vertical external forces also on those two lateral surfaces (then it all works and the uniform deformation above is correct)... but more often this is neglected.
Here we shall only discuss general relativistic diffeomorphism-invariant matter theories in a curved spacetime in the classical limit $\hbar\to 0$ for simplicity. In particular, we will not discuss the SEM pseudotensor for the gravitational field, but only the stress-energy-momentum (SEM) tensor for matter (m) fields $\Phi^A$. We emphasize that our conclusions will be independent of whether the EFEs are satisfied or not.
I) On one hand, the basic Hilbert/metric SEM tensor-density$^1$ is manifestly symmetric$^2$
$$\tag{1} \sqrt{|g|}T^{\mu\nu}~\equiv~{\cal T}^{\mu\nu}~:=~-2\frac{\delta S_m}{\delta g_{\mu\nu}},$$
cf. a comment by Lubos Motl. However, note that the basic definition (1) is not applicable to e.g. fermionic matter in a curved spacetime, cf. Section II.
Diffeomorphism invariance leads (via Noether's 2nd theorem) to an off-shell identity. Using the matter eqs. of motion (eom)
$$\tag{2} \frac{\delta S_m}{\delta \Phi^A}~\stackrel{m}{\approx}~0, $$
the corresponding Noether's 2nd identity reads$^3$
$$\tag{3} \nabla_{\mu} T^{\mu\nu}~\stackrel{m}{\approx}~0, $$
cf. e.g. Ref. 1. [Here the $\stackrel{m}{\approx}$ symbol means equality modulo matter eoms. The connection $\nabla$ is the Levi-Civita connection.] Eq. (3) serves as an important consistency check. A matter source $T^{\mu\nu}$ to the EFEs should satisfy eq. (3), cf. the (differential) Bianchi identity.
II) In the Cartan formalism, the fundamental gravitational field is not the metric tensor $g_{\mu\nu}$ but instead a vielbein $e^a{}_{\mu}$. The generalized Hilbert SEM tensor-density is defined as
$$\tag{4}{\cal T}^{\mu}{}_{\nu}~:=~{\cal T}^{\mu}{}_a~e^a{}_{\nu}, \qquad
{\cal T}^{\mu}{}_a ~:=~- \frac{\delta S_m}{\delta e^a{}_{\mu}}, $$
which is no longer manifestly symmetric, cf. e.g. Ref. 2.
Next we have two symmetries: local Lorentz symmetry and diffeomorphism invariance.
Firstly, local Lorentz symmetry leads (via Noether's 2nd theorem) to an off-shell identity. Using the matter eoms (2), the corresponding Noether's 2nd identity reads
$$\tag{5} {\cal T}^{\mu\nu}~\stackrel{m}{\approx}~{\cal T}^{\nu\mu},$$
i.e the generalized Hilbert SEM tensor-density (4) is still symmetric when the matter eoms are satisfied.
Secondly, diffeomorphism invariance leads (via Noether's 2nd theorem) to an off-shell identity$^4$
$$\tag{6} d_{\mu}{\cal T}^{\mu}{}_{\nu}
~=~{\cal T}^{\mu}{}_a ~d_{\nu}e^a{}_{\mu}
-\frac{\delta S_m}{\delta \Phi^A}~d_{\nu}\Phi^A. $$
Not surprisingly, eqs. (5), (6), and $(\nabla_{\nu} e)^a{}_{\mu}=0$ imply eq. (3).
III) On the other hand, the canonical SEM tensor-density
$$\tag{7} \Theta^{\mu}{}_{\nu}
~:=~\frac{\partial {\cal L}_m}{\partial\Phi^A_{,\mu}}\Phi^A_{,\nu} -\delta^{\mu}_{\nu}{\cal L}_m$$
is not always symmetric, cf. e.g. this Phys.SE post. The fact that the Lagrangian density ${\cal L}_m$ has no explicit spacetime dependence leads (via Noether's 1st theorem) to an off-shell identity
$$\tag{8} d_{\mu}\Theta^{\mu}{}_{\nu}
~=~-\frac{\delta S_m}{\delta e^a{}_{\mu}}~d_{\nu}e^a{}_{\mu}
-\frac{\delta S_m}{\delta \Phi^A}~d_{\nu}\Phi^A. $$
Recall that the gravitational field, the vielbein $e^a{}_{\mu}$, is not necessarily on-shell. Remember we are doing FT in curved spacetime rather than GR. As a consequence, the first term on the right-hand side of Noether's 1st identity (8) breaks the usual interpretation of Noether's 1st theorem as leading to an on-shell conservation law. [It is comforting to see that it gets restored for a constant vielbein with $d_{\nu}e^a{}_{\mu}=0$.]
IV) Eqs. (4), (6), and (8) imply that
$$\tag{9} d_{\mu}({\cal T}^{\mu}{}_{\nu}-\Theta^{\mu}{}_{\nu})~=~0.$$
One may show that there in general exists a Belinfante improvement tensor-density
$$\tag{10} {\cal B}^{\lambda\mu}{}_{\nu}~=~-{\cal B}^{\mu\lambda}{}_{\nu},$$
such that
$$\tag{11} {\cal T}^{\mu}{}_{\nu}-\Theta^{\mu}{}_{\nu}
~=~d_{\lambda}{\cal B}^{\lambda\mu}{}_{\nu},$$
cf. e.g. my Phys.SE answer here. Note that eqs. (10) and (11) are consistent with eq. (9).
V) Eq. (11) serves as an important consistency check of the Hilbert SEM tensor-density (4) vs. the canonical SEM tensor-density (7). Eq. (11) implies that the two corresponding Noether charges, the energy-momentum $4$-covectors
$$\tag{12} P_{\nu} ~:=~ \int\! d^3x~{\cal T}^0{}_{\nu} \quad\text{and}\quad \Pi_{\nu} ~:=~ \int\! d^3x~\Theta^0{}_{\nu} $$
are equal up to spatial boundary terms
$$\tag{13} P_{\nu}-\Pi_{\nu}~=~\int\! d^3x~d_i{\cal B}^{i0}{}_{\nu}~\sim~0,$$
cf. the divergence theorem.
References:
R.M. Wald, GR; Appendix E.1.
D.Z. Freedman & A. Van Proeyen, SUGRA, 2012; p. 181.
--
$^1$ A tensor-density ${\cal T}^{\mu\nu}= e T^{\mu\nu}$ is in this context just a tensor $T^{\mu\nu}$ multiplied with the density $e=\sqrt{|g|}$.
$^2$ Conventions: In this answer, we will use $(+,-,-,-)$ Minkowski sign convention. Greek indices $\mu,\nu,\lambda, \ldots,$ are so-called curved indices, while Roman indices $a,b,c, \ldots,$ are so-called flat indices.
$^3$ Note that eq. (3) is not a conservation law by itself. To get a conservation law, we need a Killing vector field, cf. e.g. my Phys.SE answer here.
$^4$ Here we have assumed that the matter fields $\Phi^A$ only carry flat or spinorial indices, cf. the setting of my Phys.SE answer here. If $\Phi^A$ also have curved indices, there will be further terms in eq. (6) proportional to the matter eoms.
Best Answer
Here is another way to think about it. It will help to make an analogy with the continuity equation for fluid flow. Suppose you have a fluid which is flowing with some non-uniform velocity $\vec{u}(\vec{r})$. Consider region $\Omega$ of the fluid bounded by a surface $S$, and assume you are told that the total mass in the region $\Omega$ is constant, and moreover the mass density $\rho(\vec{r})$ at each point in $\Omega$ is constant. What could you conclude from this? Well you know there is a mass current in the fluid given by $\vec{J}=\rho \vec{u}$, and since each piece of mass moves continuously, the change in density at a point is equal to the amount of mass flowing into that point. This is expressed mathematically by the "continuity equation" $\dot{\rho} = -\nabla \cdot \vec{J}$. Integrating this over the whole volume, we find $\dot{M} = \int_\Omega \dot{\rho} = \int_\Omega - \nabla \cdot \vec{J} = - \oint_S \vec{J} \cdot \hat{n}\, dA$. On the far right hand side the term $\vec{J} \cdot \hat{n}$ represents the mass flux through the surface.
The stress tensor is in fact not much harder to understand than this. You have two questions: Why is the force on a surface linear in $\hat{n}$? and Why is $\sigma$ symmetric? I will answer these one at a time. Both answers will be done by analogy with the fluid flow.
Why is the force on a surface linear in $\hat{n}$?
Assume we have a large piece of material, and we are told that the momentum density, which I will denote $\vec{p}$, and the angular momentum density, which I will denote $\vec{\ell}$, is constant in some region $\Omega$ with boundary $S$. Let us pick a coordinate system and pick a component of the momentum to look at, say the $i$th component, $p_i$. Then $p_i$ is analogous to $\rho$. Since each little piece of material only exerts forces on its neighbors (there are no long range forces), the $p_i$ must move continuously through the material. Thus the flow of $p_i$ is described by some current $\sigma_{ij}$ which is analogous to $-J_j$ (notice there is a sign convention). Identifying $\dot{p}_i$ with $f_i$, the force per unit volume, the continuity equation gives us $f_i = \partial_j \sigma_{ij}$ (remembering the sign convention). Integrating over the region $\Omega$, with $P_i$ being $i$th component of the total momentum, we find $\dot{P}_i = \int_\Omega \dot{p}_i = \int_\Omega \partial_j \sigma_{ij} = \oint_S \sigma_{ij} n_j \, dA$. Thus the term $\sigma_{ij} n_j$ has the interpretation of momentum flux through the surface, or in other words, the force per unit area on the surface. Thus you the answer to the first part of your question is that the force per unit area is linear in $\hat{n}$ for the same reason that mass flux through an area is linear in $\hat{n}$, and this reason is that there is some current describing how the mass (or momentum) is moving through the medium, and the flux is just the current dotted into $\hat{n}$.
Why is $\sigma$ symmetric?
Now lets address the second part of your question, why must $\sigma_{ij}$ be symmetric if the object is in equilibrium. Let us now consider the $i$th component of the angular momentum density $\ell_i$. We know that an external agent exerting a force per unit area $\vec{f}$ at a point $\vec{r}$ on the surface is exerting a torque whose $i$th component is given by $\tau_i = \epsilon_{ijk} r_j f_k$. However we know from the previous paragraph that $f_k = \sigma_{kh} n_h$. Thus we conclude that $\tau_i$, which is the flux of $\ell_i$, is given by $\epsilon_{ijk}r_j \sigma_{kh} n_h$. We know that this ought to be a current dotted into $\hat{n}$, so the $\ell_i$ current must be $\epsilon_{ijk}r_j \sigma_{kh}$. The change in $\ell_i$ per unit volume, which is $i$th component of the torque per unit volume $\tau_i$, is the divergence of this current: $\dot{\ell}_i = \partial_h \epsilon_{ijk}r_j \sigma_{kh} = \epsilon_{ijk} (\partial_h r_j) \sigma_{kh} + \epsilon_{ijk}r_j (\partial_h \sigma_{kh}) = \epsilon_{ijk} \delta_{hj} \sigma_{kh} + \epsilon_{ijk}r_j f_k = \epsilon_{ijk} \sigma_{kj} + \epsilon_{ijk}r_j f_k.$
The second term is $\vec{r} \times \vec{f}$ as expected, this takes into account the angular momentum produced by uniform translation of the small piece of material. The other term is the antisymmetric part of $\sigma$ and it represents a rotation of the small piece of material about its center of mass. To show that $\sigma$ must be symmetric at an arbitrary point $\vec{r}$ we first move the origin to $\vec{r}$ and then find the expression for $\dot{\ell}_i$, which must be zreo since the object is in equillibrium. We find $0 = \dot{\ell}_i = \epsilon_{ijk} \sigma_{kj}$ where the $\vec{r} \times \vec{f}$ term was dropped because $\vec{r}$ is zero. Thus we find that the antisymmetric part of $\sigma$ must be 0.