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.
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
I will answer your questions below.
Yes. The Cauchy stress is still symmetric. The thing is that this stress measure is not energy conjugate with the strain measure that you picked. This means that the double contraction between your stress and the strain (rate) tensor does not represent the work (rate) done to achieve that deformation level.
It is still symmetric but you should pick another pair of strain-stress. The "right" stress tensor for your strain measure is the second Piola-Kirchoff stress tensor. I think that the Cauchy stress tensor does not have a conjugate pair, though (see 1).
That depends on the material that you are using. My first guess is that the football can be modeled with a hyperelastic material (see 2).
References
Hoger, A. (1987). The stress conjugate to logarithmic strain. International Journal of Solids and Structures, 23(12), 1645-1656.
Karimi, A., Razaghi, R., Navidbakhsh, M., Sera, T., & Kudo, S. (2016). Measurement of the mechanical properties of soccer balls using digital image correlation method. Sport Sciences for Health, 12(1), 69-76.