Where is the Elastic Strain Energy per Actual Volume Gone? I have a probably real stupid question, but I still can not wrap my head around it. 
The change in inner energy per unit volume $v$ of the actual state of a stressed crystal is given by the well known formula
$$
de_v=\sigma_{ij}d\epsilon_{ij} \label{eq1} \tag{1}
$$
If the crystal is only homogeneously deformed, this formula should simplify to
$$
de_v=\sigma d\epsilon \label{eq2} \tag{2}
$$
With $\sigma=\sigma_{kk}$ and $\epsilon=\epsilon_{kk}$ (using Einsteins summation convention).
I know how to derive this formula using virtual displacements, but I wondered if I can also derive this formula from the general formula for internal energy with volume work performed by the system.
So, I start with the mother of all thermodynamical equations for a single phase system
$$
e=-vp+TS+\mu N \label{eq3} \tag{3}
$$
$v$ being the actual volume and $p$ the hydrostatical pressure, $e$ the internal energy, $T$ temperature, $S$ entropy, $\mu$ chemical potential and $N$ the number of particles. The Gibbs-Duhem Relation for a process reads
$$
-vdp+SdT+Nd\mu=0 \label{eq4} \tag{4}
$$
After deviding \eqref{eq4} by the actual volume of the system $v$, we end up with
$$
-dp+s_v dT+\rho d\mu = 0 \label{eq5} \tag{5}
$$
Here the subscript $v$ means the value per volume and $\rho = N/v$.
I also devide equation \eqref{eq3} by the systems actual volume $v$ and calculate the total differential of the resulting equation.
$$
de_v = -dp + s_v dT + Tds_v+\rho d\mu + \mu d\rho \label{eq6} \tag{6}
$$
Using \eqref{eq5}, I immediately obtain
$$
de_v =  T ds_v + \mu d\rho \label{eq7} \tag{7}
$$
So, the internal energy density seems to be independent of the deformation of the crystal. There is no term that could be understood as resembling the $de_v=\sigma d\epsilon$ term from equation \eqref{eq2}.
What am I missing????!!!! 
 A: The reason that you analysis is not correct is that your equation for the internal energy per unit volume is not correct.  The internal energy per unit volume is a function not only of the deformation but also of one other intensive variable.  The reference I provided very nicely lays this all out for a Hookean (linearly elastic) solid.
A: OK.  Here's my spin on this.  I am going to assume, as in the other analyses this thread, that the loading is one of isotropic stress and that we are dealing with small volumetric strains.  So the basic equation for change in specific internal energy is given by:
$$du=Tds+v_0\sigma d\epsilon\tag{1}$$where $\epsilon$ is the volumetric strain, $\sigma$ is the isotropic stress, and $v_0$ is the (reference) specific volume when $\sigma = 0$, $\epsilon = 0$, and $T = T_0$.  In addition to Eqn. 1, the other key equation needed is the equation of state (EOS)for the ideal Hookean elastic solid:
$$\epsilon=\alpha (T-T_0)+\frac{\sigma}{K}\tag{2}$$where $\alpha$ is the coefficient of volumetric thermal expansion, and K is the bulk modulus of the solid:  $$K=\frac{E}{3(1-2\nu)}$$where E is the Young's modulus and $\nu$ is the Poisson's ratio.
From Eqn. 1 and the definition of the Helmholtz free energy a, we have for the differential of a:
$$da=-sdT+v_0\sigma d\epsilon\tag{3}$$
In Eqn. 1, if we substitute $$ds=\left(\frac{\partial s}{\partial T}\right)_{\epsilon}dT+\left(\frac{\partial s}{\partial \epsilon}\right)_{T}d\epsilon$$we obtain:
$$du=C_vdT+\left[v_0\sigma+T\left(\frac{\partial s}{\partial \epsilon}\right)_{T}\right]d\epsilon\tag{5}$$where $$C_v=T\left(\frac{\partial s}{\partial T}\right)_{\epsilon}$$
From Eqn. 3, we can derive the following Maxwell equation for the partial derivative of the specific entropy s with respect to the strain:
$$\left(\frac{\partial s}{\partial \epsilon}\right)_{T}=-v_0\left(\frac{\partial \sigma}{\partial T}\right)_{\epsilon}\tag{6}$$Substituting the EOS Eqn. 2 into Eqn. 6 then yields:  $$\left(\frac{\partial s}{\partial \epsilon}\right)_{T}=v_0K\alpha\tag{7}$$
Substitution of Eqn. 7 into Eqn. 5 then yields:
$$du=C_vdT+v_0(\sigma+\alpha KT)d\epsilon\tag{8}$$
If we next substitute EOS Eqn. 2 for $epsilon$ into Eqn. 8, we obtain:
$$du=C_vdT+v_0(K\alpha^2TdT+\frac{\sigma d\sigma}{K}+\alpha d(\sigma T))\tag{9}$$Eqn. 9 is an exact differential for du, and integrates immediately to:
$$u=u_0+C_v(T-T_0)+v_0\left[K\alpha^2\frac{(T^2-T_0^2)}{2}+\frac{\sigma^2}{2K}+\alpha T\sigma\right]\tag{10}$$
If we apply similar procedures to the entropy s, we obtain:
$$s=s_0+C_v\ln{(T/T_0)}+v_0\left[K\alpha^2(T-T_0)+\alpha T\right]\tag{11}$$
Eqns. 10 and 11 can be combined to obtain the Helmholtz free energy (to quadratic terms in $(T-T_0)$ as follows:  $$a=u-Ts=a_0-s_0(T-T_0)-\frac{C_v(T-T_0)^2}{2T_0}-\frac{v_0K\alpha^2}{2}(T-T_0)^2+\frac{v_0\sigma^2}{2K}\tag{12}$$
This analysis can be continued one more step to provide an explicit relationship for the specific Gibbs free energy g (aka the chemical potential $\mu$) by subtracting $v_0\sigma \epsilon $, and illustrating how the chemical potential is related to the stress $\sigma$ and the temperature T.
A: Drawing from a previous comment from @Chemomechanics, I want to post one possible way out of this dilemma. 
First of all, the $\rho$ in equation (7) (from the question) is, of course, not independent of the volume but it holds that $\rho = N/V$. We can thus calculate the total differential
$$
d\rho = \frac{1}{V} dN - \frac{N}{V^2}dV \tag{8} \label{eq8}
$$
Inserting into equation (7)(from the question) results in
$$
de_v = Tds_v + \frac{\mu}{V} dN - \frac{\mu N}{V^2}dV \tag{9} \label{eq9}
$$
To see how this relates to $de_v=\sigma d\epsilon$ we'll have to calculate the chemical potential $\mu$. Lets first, for the sake of simplicity assume that $N$ is constant and that the process is isothermal, that is $dT=0$. From the Gibbs-Duhem equation (4)(from the question) it follows, that
$$
d\mu = \frac{V}{N}dp \tag{10} \label{eq10}
$$
For small strains we can now assume a constitutive relation of the form
$$
p(V)=-K(\frac{V}{V_0}-1)\tag{11} \label{eq11}
$$
$K$ being the Bulk-Modulus and $V_0$ the volume at which $p=0$. Inverting \eqref{eq11} and integrating \eqref{eq10} from $0$ to $p$, leads to
$$
\mu(N,p)=\frac{p}{2N}(V+V_0) \tag{12} \label{eq12}
$$
Inserting this, $N=const$ and $dV/V=d\epsilon$ into \eqref{eq9}, we end up with
$$
de_v = Tds_v - \frac{1}{2} p ( 1 + \frac{V_0}{V}) d\epsilon \tag{13} \label{eq13}
$$
Since in the small strain approximation $V_0 \approx V$, we got
$$
de_v = Tds_v - p d\epsilon \tag{14} \label{eq14}
$$
Identifying $-p$ with $\sigma$, we are left with
$$
de_v = Tds_v + \sigma d\epsilon \tag{15} \label{eq15}
$$
This way, we can recover the dependence of the internal energy per volume on the strain.
A: I think I now have found a more general explanation then the previously given one.  ( I was not sure weather to make it an update or a new answer, but since it is quite long, I made it a new one. I can change that if it needs to be.. )
I would like to start with the observation that this problem actualy really is a problem and did not emerge because I did forget any variables ( or dependencies ). I actually saw dependencies where there were none. But slowly.
We again start with the first law of thermodynamics, writing
$$
de=-pdv+\mu dN + Tds \tag{16} \label{eq16}
$$
This is the differential for the function $e=e(v,N,s)$. Because $e$ is themordynamical potential, the Euler equation 
$$
e(\lambda v,\lambda N,\lambda s)=\lambda e(v,N,s) \tag{17} \label{eq17}
$$
holds. But since $N=\rho v$ and $S=s_v v$, we can eliminate the dependence of $e$ on $v$ by using \eqref{eq17}
$$
e(v,v \rho,s_v v) = v e(1,\rho,s)=:v e_v (\rho, s_v) \tag{18} \label{eq18}
$$
I am doing this seemingly useless calculation to explicitly show, that $e_v$ does not depend on the volume $v$. With the densities $\rho$ and $s_v$, we introduced new natural variables for the function $e_v=e_v(\rho,s_v)$. 
So we are really left with the following differential for $e_v$ ( following as in equation \eqref{eq7} )
$$
de_v= Tds_v + \mu d\rho \tag{19} \label{eq19}
$$
So all the information on the volume, and hence on the strain, is incorporated into those various densities. To extract explicitly the dependence on the strain we actually have to make use of the small strain approximation.
We introduce a constant reference volume $v_r$, for which we can assume that the small strain approximation holds. That is, we postulate that there exists a volume $v_r$ so that for all volumes $v$ the system will be exhibiting, the deformation relative to $v_r$ is small. That is
$$
\frac{v_r}{v}=\frac{v+dv}{v}=(1+\epsilon) =:J \approx 1 \tag{20} \label{eq20}
$$
We then define the densities $\rho_r = N / v_r = \rho / J$, $s_r = S / v_r = s_v / J$ and $e_r = e / v_r = e_v / J$. We furthermore calculate the various differentials that we gonna need using the small strain approximation.
$$
dJ = - \frac{v_r}{v^2}dV = - J d\epsilon \tag{21} \label{eq21}
$$
$$
d \rho = d(\rho_r J) = J d\rho_r - \rho_r J d\epsilon = d\rho_r - \rho_r d\epsilon \tag{22} \label{eq22}
$$
$$
d s_v = d(s_r J) = J ds_r - s_r J d\epsilon = ds_r - s_r d\epsilon \tag{23} \label{eq23}
$$
$$
d e_v = d(e_r J) = J de_r - e_r J d\epsilon = de_r - e_r d\epsilon \tag{24} \label{eq24}
$$
Inserting \eqref{eq22} to \eqref{eq24} into \eqref{eq19}, we end up with
$$
de_r - e_r d\epsilon = Tds_r - Ts_r d\epsilon + \mu d\rho_r - \mu \rho_r d\epsilon 
$$
$$
de_r = Tds_r + \mu d \rho_r + (e_r-Ts_r - \mu \rho_r)
d \epsilon \tag{25} \label{eq25}
$$
Deviding the equation for $e$, \eqref{eq3}, by $v_r$ and using the small strain approximation again, we find, that
$$
e_r = T s_r + \mu \rho_r - p J \approx T s_r + \mu \rho_r - p
$$
Or 
$$
p = e_r - T s_r - \mu \rho_r \tag{26} \label{eq26}
$$
Using this in equation \eqref{eq25}, we finally arrive at
$$
de_r = T ds_r+ \mu d \rho_r - p d\epsilon = Tds_r + \mu d \rho_r + \sigma d\epsilon \tag{27} \label{eq27}
$$
So, to regain the dependencie of $e_v$ on $\epsilon$, we have to use densities that were calculated relative to a constant reference volume ( which we could have done from the bloody beginning ), that is $e_r$, $s_r$ and $\rho_r$. If we use $e_v$, any changes in strain have to be incorporated into the densities.
