应变能密度的推导
考虑求解域 \(V\) ，弹性力学问题的平衡方程与边界条件为 \[ \tag{1} \begin{cases} \nabla\cdot\boldsymbol{\sigma}+\boldsymbol{b}=\rho\ddot{\boldsymbol{u}}&~\mathrm{in}~V\\[5pt] \boldsymbol{n}\cdot\boldsymbol{\sigma}=\boldsymbol{t}&~\mathrm{on}~\partial V\\[5pt] \boldsymbol{u}=\boldsymbol{\bar{u}}&~\mathrm{on}~\partial V^u \end{cases} \]
式中，\(\rho\) 为密度，\(\boldsymbol{b}\) 为体力，\(\boldsymbol{t}\) 为自然边界 \(\partial V\) 上的面力，\(\boldsymbol{\bar{u}}\) 为强制边界 \(\partial V^u\) 上的位移。
由热力学第一定律可知 \[ \tag{2} \delta V_\varepsilon+\delta K=\delta W \] 式中，\(\delta V_\varepsilon\)、\(\delta K\)、 \(\delta W\) 分别为应变能（势能）增量、动能增量、外力功增量。
应变能增量
弹性应变能与变形过程无关，仅取决于变形的最终状态，即为状态函数，其增量为全微分 \[ \tag{3a} v_\varepsilon=v_\varepsilon(\boldsymbol{\varepsilon}) \]
\[ \tag{3b} \delta v_\varepsilon=\frac{\partial v_\varepsilon}{\partial \boldsymbol{\varepsilon}}:\delta\boldsymbol{\varepsilon} \]
式中 \(v_\varepsilon\) 为应变能密度。于是可得到应变能增量 \[ \tag{4} \delta V_\varepsilon=\int_V\delta v_\varepsilon\mathrm{d}v=\int_V\frac{\partial v_\varepsilon}{\partial \boldsymbol{\varepsilon}}:\delta\boldsymbol{\varepsilon}\mathrm{d}v \]
动能增量
动能 \(K\) 等于 \[ \tag{5} K(\dot{\boldsymbol{u}})=\int_V\frac{1}{2}\rho\dot{\boldsymbol{u}}\cdot\dot{\boldsymbol{u}}\mathrm{d}v \] 于是动能增量为 \[ \tag{6} \begin{aligned} \delta K&=K(\dot{\boldsymbol{u}}+\delta\dot{\boldsymbol{u}})-K(\dot{\boldsymbol{u}})\\[5pt] &=\int_V\frac{1}{2}\rho(\dot{\boldsymbol{u}}+\delta\dot{\boldsymbol{u}})\cdot(\dot{\boldsymbol{u}}+\delta\dot{\boldsymbol{u}})\mathrm{d}v-\int_V\frac{1}{2}\rho\dot{\boldsymbol{u}}\cdot\dot{\boldsymbol{u}}\mathrm{d}v\\[5pt] &=\int_V\big[\rho\dot{\boldsymbol{u}}\cdot\delta\dot{\boldsymbol{u}}+\underbrace{\bcancel{o(\delta^2\dot{\boldsymbol{u}})}}_{\text{高阶无穷小}}\big]\mathrm{d}v=\int_V\rho(\dot{\boldsymbol{u}}\cdot\ddot{\boldsymbol{u}})\delta t\mathrm{d}v\\[1pt] &=\int_V\rho\ddot{\boldsymbol{u}}\cdot\delta\boldsymbol{u}\mathrm{d}v \end{aligned} \]
外力功增量
对于任一微元体 \(V\)，其表面为 \(\partial V\)，外载荷在物体上做功为 \[ \tag{7} \begin{aligned} \delta W&=\int_V \boldsymbol{b}\cdot\delta\boldsymbol{u}\mathrm{d}v+\oint_{\partial V}\boldsymbol{t}\cdot\delta\boldsymbol{u}\mathrm{d}s\\[8pt] &=\int_V \boldsymbol{b}\cdot\delta\boldsymbol{u}\mathrm{d}v+\oint_{\partial V}\boldsymbol{n}\cdot\boldsymbol{\sigma}\cdot\delta\boldsymbol{u}\mathrm{d}s \end{aligned} \]
由高斯公式
\[ \oint_{\partial V}\boldsymbol{n}\cdot\boldsymbol{F}\mathrm{d}s=\int_V\nabla\cdot\boldsymbol{F}\mathrm{d}v \]
式中，向量场 \(\boldsymbol{F}=P\boldsymbol{i}+Q\boldsymbol{j}+R\boldsymbol{k}\)
可得
\[ \tag{8} \delta W=\int_V \boldsymbol{b}\cdot\delta\boldsymbol{u}\mathrm{d}v+\int_V\nabla\cdot(\boldsymbol{\sigma}\cdot\delta\boldsymbol{u})\mathrm{d}v \] 考虑到 \[ \tag{9} \begin{aligned} \nabla\cdot(\boldsymbol{\sigma}\cdot\delta\boldsymbol{u})&=\partial_i\boldsymbol{e}_i\cdot(\sigma_{jk}\boldsymbol{e}_j\boldsymbol{e}_k\cdot\delta u_l\boldsymbol{e}_l)=\delta_{ij}\delta_{kl}(\sigma_{jk}\delta u_l)_{,i}\\[5pt] &=(\sigma_{il}\delta u_l)_{,i}=\sigma_{il,i}\delta u_l+\sigma_{il}\delta u_{l,i} \end{aligned} \]
进一步
\[ \tag{10} \begin{aligned} \sigma_{il}\delta u_{l,i}&=\frac{1}{2}\sigma_{il}\delta u_{l,i}+\frac{1}{2}\sigma_{li}\delta u_{i,l}\\[5pt] &=\frac{1}{2}\sigma_{il}\delta u_{l,i}+\frac{1}{2}\sigma_{il}\delta u_{i,l}\\[5pt] &=\sigma_{il}\frac{1}{2}(\delta u_{l,i}+\delta u_{i,l})=\sigma_{il}\delta \varepsilon_{il} \end{aligned} \]
式\((10)\)的推导利用了应力张量 \(\boldsymbol{\sigma}\) 的对称性和应变张量 \(\boldsymbol{\varepsilon}\) 的定义，即 \(\boldsymbol{\varepsilon}=[\nabla\otimes\boldsymbol{u}+\boldsymbol{u}\otimes\nabla]/2\) 。于是有 \[ \tag{11} \nabla\cdot(\boldsymbol{\sigma}\cdot\delta\boldsymbol{u})=(\nabla\cdot\boldsymbol{\sigma})\cdot\delta\boldsymbol{u}+\boldsymbol{\sigma}:\delta \boldsymbol{\varepsilon} \] 将式\((11)\)代入式\((8)\)可得 \[ \tag{12} \delta W=\int_V(\nabla\cdot\boldsymbol{\sigma}+\boldsymbol{b})\cdot\delta\boldsymbol{u}\mathrm{d}v+\int_V\boldsymbol{\sigma}:\delta\boldsymbol{\varepsilon}\mathrm{d}v \]
将式\((4)\)、\((6)\)、\((12)\)代入式\((2)\)可得 \[ \tag{13} \delta W-\delta V_\varepsilon-\delta K=\int_V(\nabla\cdot\boldsymbol{\sigma}+\boldsymbol{b}-\rho\ddot{\boldsymbol{u}})\cdot\delta\boldsymbol{u}\mathrm{d}v+\int_V\big(\boldsymbol{\sigma}-\frac{\partial v_\varepsilon}{\partial \boldsymbol{\varepsilon}}\big):\delta\boldsymbol{\varepsilon}\mathrm{d}v=0 \] > \(\delta\boldsymbol{u}\) 与 \(\delta\boldsymbol{\varepsilon}\) 并不独立，也就无法由式\((13)\)导出式\((1)\)中的平衡方程，以及式\((14)\)中的物理方程 > > 式\((1)\)中的平衡方程可由牛顿第二定律证明
由式\((1)\)中的平衡方程可知，式\((13)\)中 \(\nabla\cdot\boldsymbol{\sigma}+\boldsymbol{b}-\rho\ddot{\boldsymbol{u}}=\boldsymbol{0}\) ，于是有 \[ \tag{14a} \int_V\big(\boldsymbol{\sigma}-\frac{\partial v_\varepsilon}{\partial \boldsymbol{\varepsilon}}\big):\delta\boldsymbol{\varepsilon}\mathrm{d}v=0 \] 即有 \[ \tag{14b} \boldsymbol{\sigma}=\frac{\partial v_\varepsilon}{\partial \boldsymbol{\varepsilon}} \]
式\((14\mathrm{b})\)被称为Green公式，即能量形式的物理方程。推导式\((14)\)的过程中使用了热力学第一定律、高斯公式、应力张量的对称性、应变张量的定义、平衡方程，适用于一般的弹性材料，不局限线弹性材料。
对于线弹性材料，积分与路径无关，可假设按等比例加载，则应变能密度为 \[ \tag{15} \begin{aligned} v_\varepsilon&=\int_{t_0}^{t_1}\boldsymbol{\sigma}:\delta\boldsymbol{\varepsilon}=\int_0^1t\boldsymbol{\sigma}:(t\delta\boldsymbol{\varepsilon})\\[5pt] &=\boldsymbol{\sigma}:\boldsymbol{\varepsilon}\int_0^1t\delta t=\frac{1}{2}\boldsymbol{\sigma}:\boldsymbol{\varepsilon}=\frac{1}{2}\sigma_{ij}\varepsilon_{ij} \end{aligned} \] 对于各向同性线弹性材料，其 \(\mathrm{Lam} \acute {\mathrm{e}}\) 形式的物理方程为 \[ \tag{16a} \boldsymbol{\sigma}(\boldsymbol{\varepsilon})=2\mu\boldsymbol{\varepsilon}+\lambda\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{I} \] 其应变能密度为 \[ \tag{16b} v_\varepsilon(\boldsymbol{\varepsilon})=\frac{1}{2}\boldsymbol{\sigma}:\boldsymbol{\varepsilon}=\mu\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}+\frac{1}{2}\lambda\mathrm{tr}^2(\boldsymbol{\varepsilon})=\mu\mathrm{tr}(\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon})+\frac{1}{2}\lambda\mathrm{tr}^2(\boldsymbol{\varepsilon}) \]
应变能密度还可由应变张量的三个主应变（特征值）给出 \[ \tag{16c} v_\varepsilon(\varepsilon_1,\varepsilon_2,\varepsilon_3)=\mu(\varepsilon_1^2+\varepsilon_2^2+\varepsilon_3^2)+\frac{1}{2}\lambda(\varepsilon_1+\varepsilon_2+\varepsilon_3)^2 \]
\(\mathrm{Hooke}\) 形式的物理方程
\[ \boldsymbol{\varepsilon}(\boldsymbol{\sigma})=\frac{1+\nu}{E}\boldsymbol{\sigma}-\frac{\nu}{E}\mathrm{tr}(\boldsymbol{\sigma})\boldsymbol{I} \] 应变能密度 \[ v_\varepsilon(\boldsymbol{\sigma})=\frac{1}{2}\boldsymbol{\sigma}:\boldsymbol{\varepsilon}=\frac{1+\nu}{2E}\boldsymbol{\sigma}:\boldsymbol{\sigma}-\frac{\nu}{2E}\mathrm{tr}^2(\boldsymbol{\sigma}) \] 主应力形式 \[ v_\varepsilon(\sigma_1,\sigma_2,\sigma_3)=\frac{1+\nu}{2E}(\sigma_1^2+\sigma_2^2+\sigma_3^2)-\frac{\nu}{2E}(\sigma_1+\sigma_2+\sigma_3)^2 \]
应变能密度的分解
首先将应变张量 \(\boldsymbol{\varepsilon}\) 分解为球形部分和偏斜部分 \[ \tag{17} \begin{cases} \boldsymbol{\varepsilon}_\mathrm{v}=\displaystyle\frac{1}{3}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{I}\\[5pt] \boldsymbol{\varepsilon}_\mathrm{d}=\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_\mathrm{v} \end{cases} \] 不难证明 \(\boldsymbol{\varepsilon}_\mathrm{v}\) 和 \(\boldsymbol{\varepsilon}_\mathrm{d}\) 满足下列表达式 \[ \tag{18} \begin{cases} \boldsymbol{\varepsilon}_\mathrm{v}:\boldsymbol{\varepsilon}_\mathrm{d}=0\\[5pt] \boldsymbol{\varepsilon}_\mathrm{v}:\boldsymbol{\varepsilon}_\mathrm{v}=\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}_\mathrm{v}=\displaystyle\frac{1}{3}\mathrm{tr}^2(\boldsymbol{\varepsilon})\\[5pt] \boldsymbol{\varepsilon}_\mathrm{d}:\boldsymbol{\varepsilon}_\mathrm{d}=\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}_\mathrm{d}=\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}-\displaystyle\frac{1}{3}\mathrm{tr}^2(\boldsymbol{\varepsilon}) \end{cases} \] 考虑到 \(\boldsymbol{\varepsilon}\)、\(\boldsymbol{\varepsilon}_\mathrm{v}\)、\(\boldsymbol{\varepsilon}_\mathrm{d}\) 三者的对称性，有 \[ \tag{19} \begin{cases} \boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}_\mathrm{v}=\displaystyle\frac{1}{3}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{\varepsilon}\\[5pt] \boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}_\mathrm{d}=\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}-\displaystyle\frac{1}{3}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{\varepsilon}\\[5pt] \boldsymbol{\varepsilon}_\mathrm{v}\cdot\boldsymbol{\varepsilon}_\mathrm{d}=\displaystyle\frac{1}{3}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{\varepsilon}-\frac{1}{9}\mathrm{tr}^2(\boldsymbol{\varepsilon})\boldsymbol{I}\\[8pt] \boldsymbol{\varepsilon}_\mathrm{v}\cdot\boldsymbol{\varepsilon}_\mathrm{v}=\displaystyle\frac{1}{9}\mathrm{tr}^2(\boldsymbol{\varepsilon})\boldsymbol{I}\\[5pt] \boldsymbol{\varepsilon}_\mathrm{d}\cdot\boldsymbol{\varepsilon}_\mathrm{d}=\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}-\displaystyle\frac{2}{3}\mathrm{tr}(\boldsymbol{\varepsilon})\boldsymbol{\varepsilon}+\frac{1}{9}\mathrm{tr}^2(\boldsymbol{\varepsilon})\boldsymbol{I} \end{cases} \] 将应变张量的分解形式代入式\((16\mathrm{a})\)可得 \[ \tag{20} \begin{aligned} \boldsymbol{\sigma}&=2\mu(\boldsymbol{\varepsilon}_\mathrm{v}+\boldsymbol{\varepsilon}_\mathrm{d})+\lambda\mathrm{tr}(\boldsymbol{\varepsilon}_\mathrm{v}+\boldsymbol{\varepsilon}_\mathrm{d})\boldsymbol{I}\\[5pt] &=\underbrace{(2\mu+3\lambda)\boldsymbol{\varepsilon}_\mathrm{v}}_{\boldsymbol{\sigma}_\mathrm{v}}+\underbrace{2\mu\boldsymbol{\varepsilon}_\mathrm{d}}_{\boldsymbol{\sigma}_\mathrm{d}} \end{aligned} \] 式\((20)\)定义了应力张量的球形部分和偏斜部分，且不难证明 \(\boldsymbol{\varepsilon}_\mathrm{v}\xrightarrow{\text{式(16a)}}\boldsymbol{\sigma}_\mathrm{v}\) 、\(\boldsymbol{\varepsilon}_\mathrm{d}\xrightarrow{\text{式(16a)}}\boldsymbol{\sigma}_\mathrm{d}\) 。
将式\((20)\)代入式\((15)\)，并利用式\((18)\)可得 \[ \tag{21} \begin{aligned} v_\varepsilon&=\frac{1}{2}[(2\mu+3\lambda)\boldsymbol{\varepsilon}_\mathrm{v}+2\mu\boldsymbol{\varepsilon}_\mathrm{d}]:(\boldsymbol{\varepsilon}_\mathrm{v}+\boldsymbol{\varepsilon}_\mathrm{d})\\[8pt] &=\underbrace{\frac{1}{6}(2\mu+3\lambda)\mathrm{tr}^2(\boldsymbol{\varepsilon})}_{v_\mathrm{v}}+\underbrace{\mu\big[\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}-\frac{1}{3}\mathrm{tr}^2(\boldsymbol{\varepsilon})\big]}_{v_\mathrm{d}} \end{aligned} \] 式\((21)\)定义了应变能密度的球形部分和偏斜部分，且不难证明 \(\boldsymbol{\varepsilon}_\mathrm{v}\xrightarrow{\text{式(16b)}}v_\mathrm{v}\) 、\(\boldsymbol{\varepsilon}_\mathrm{d}\xrightarrow{\text{式(16b)}}v_\mathrm{d}\) 。
一般将 \(v_\mathrm{v}\) 称为体积改变能密度，\(v_\mathrm{d}\) 称为畸变能密度。
采用应变张量的三个主应变， \(v_\mathrm{v}\)、\(v_\mathrm{d}\) 可写作 \[ \tag{22} \begin{cases} v_\mathrm{v}=\displaystyle\frac{1}{6}(2\mu+3\lambda)(\varepsilon_1+\varepsilon_2+\varepsilon_3)^2\\[8pt] v_\mathrm{d}=\displaystyle\frac{1}{3}\mu\big[(\varepsilon_1-\varepsilon_2)^2+(\varepsilon_2-\varepsilon_3)^2+(\varepsilon_3-\varepsilon_1)^2\big] \end{cases} \]
若将应力张量 \(\boldsymbol{\sigma}\) 分解为球形部分和偏斜部分 \[ \begin{cases} \boldsymbol{\sigma}_\mathrm{v}=\displaystyle\frac{1}{3}\mathrm{tr}(\boldsymbol{\sigma})\boldsymbol{I}\\[5pt] \boldsymbol{\sigma}_\mathrm{d}=\boldsymbol{\sigma}-\boldsymbol{\sigma}_\mathrm{v} \end{cases} \] 采用 \(\mathrm{Hooke}\) 形式的物理方程，则相应的结论为 \[ \boldsymbol{\varepsilon}=\underbrace{\frac{1-2\nu}{E}\boldsymbol{\sigma}_\mathrm{v}}_{\boldsymbol{\varepsilon}_\mathrm{v}}+\underbrace{\frac{1+\nu}{E}\boldsymbol{\sigma}_\mathrm{d}}_{\boldsymbol{\varepsilon}_\mathrm{d}} \] 应变能密度 \[ v_\varepsilon=\underbrace{\frac{1-2\nu}{6E}\mathrm{tr}^2(\boldsymbol{\sigma})}_{v_\mathrm{v}}+\underbrace{\frac{1+\nu}{2E}\big[\boldsymbol{\sigma}:\boldsymbol{\sigma}-\frac{1}{3}\mathrm{tr}^2(\boldsymbol{\sigma})\big]}_{v_\mathrm{d}} \] 主应力形式 \[ \begin{cases} v_\mathrm{v}=\displaystyle\frac{1-2\nu}{6E}(\sigma_1+\sigma_2+\sigma_3)^2\\[8pt] v_\mathrm{d}=\displaystyle\frac{1+\nu}{6E}\big[(\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2\big] \end{cases} \]