邱小伙的博客

在给出 Hashin-Shtrikman 变分原理的陈述之后，本文档关注如何构造一个极化应力场，并得到等效模量上下界的估计。

带约束的 Hashin-Shtrikman 变分原理

H-S 给出的泛函形式如下所示： \[ \begin{equation}\label{eq:hs_vari} \mathcal{H}(\boldsymbol{p},\tilde{\boldsymbol{\varepsilon}}\mid\mathbb{L}_{0}) = \int_{\Omega} \left( \frac{1}{2}\boldsymbol{p}:\tilde{\boldsymbol{\varepsilon}} + \boldsymbol{p}:\boldsymbol{\varepsilon}_{0} - \frac{1}{2} \boldsymbol{p}:(\mathbb{L} - \mathbb{L}_{0})^{-1}:\boldsymbol{p} \right) \mathrm{d}\Omega + \frac{1}{2} \boldsymbol{\varepsilon}_{0} : \mathbb{L}_{0} : \boldsymbol{\varepsilon}_{0}. \end{equation} \] 式中，泛函 \(\eqref{eq:hs_vari}\) 含有两个依赖变量：极化应力 \(\boldsymbol{p}\) 和应变场扰动 \(\tilde{\boldsymbol{\varepsilon}}\)，以及一个参数：参考模量 \(\mathbb{L}_{0}\)。依赖变量之间并不独立，而是通过如下辅助方程相互联系： \[ \begin{equation} \left( \mathbb{L}_{0} : \nabla_{s}\tilde{\boldsymbol{u}} + \boldsymbol{p} \right)\cdot\nabla = \boldsymbol{0}, \quad \tilde{\boldsymbol{u}}|_{\partial\Omega} = \boldsymbol{0}. \label{eq:subsd} \end{equation} \]

极化应力的可行域是可积函数空间；应变场扰动的可行域空间是所有相容且在区域 \(\Omega\) 内体积平均等于零的应变场；宏观应变 \(\boldsymbol{\varepsilon}_{0}\) 是给定值。泛函 \(\eqref{eq:hs_vari}\) 在如下位置取驻值： \[ \begin{equation} \boldsymbol{p} = \Delta \mathbb{L}:(\boldsymbol{\varepsilon}_{0} + \tilde{\boldsymbol{\varepsilon}}). \label{eq:stat} \end{equation} \] 当参考模量 \(\mathbb{L}_{0}^{-}\) 使得 \(\Delta\mathbb{L} \triangleq \mathbb{L} - \mathbb{L}_{0}\) 正定时，泛函 \(\eqref{eq:hs_vari}\) 在驻点处取最大值；而当参考模量 \(\mathbb{L}_{0}^{+}\) 使 \(\Delta\mathbb{L}\) 负定时，泛函在驻点处取最小值。因此有如下不等式成立： \[ \mathcal{H}(\hat{\boldsymbol{p}}^{-},\tilde{\boldsymbol{\varepsilon}}\mid\mathbb{L}_{0}^{-}) \leq \frac{1}{2} \boldsymbol{\varepsilon}_{0} : \mathbb{L}^{c} : \boldsymbol{\varepsilon}_{0} \leq \mathcal{H}(\hat{\boldsymbol{p}}^{+},\tilde{\boldsymbol{\varepsilon}}\mid\mathbb{L}_{0}^{+}), \] 式中，\(\hat{\boldsymbol{p}}\) 是任意满足可积条件的极化应力场。接下来将构造分片常值分布的极化应力场，并显式地写出泛函取值。

极化应力场的构造

以下将给出 Hashin-Shtrikman 变分原理的表述。构造用于该变分原理的试验场将在其它文档中给出。

双场变分原理的表述

以下将通过 Young-Fechel 变换获得 H-S 变分原理中的双场泛函形式。考虑区域 \(\Omega\) 内的非均质材料，模量为 \(\mathbb{L}\)。给定位移边界条件 \(\boldsymbol{u}|_{\partial\Omega}=\boldsymbol{u}_{0}\)，由此得到在该区域内的最小势能原理表述为 \[ \check\Pi = \min_{\boldsymbol{u}\in \mathcal{U}} \Pi(\boldsymbol{u}) = \min_{\boldsymbol{u}\in \mathcal{U}}\frac{1}{2}\int_{\Omega} \nabla_{s}\boldsymbol{u}:\mathbb{L}:\nabla_{s}\boldsymbol{u} \ \mathrm{d}\Omega, \quad \mathcal{U} \triangleq \{ \boldsymbol{u} \mid \boldsymbol{u} \in H^{1},\ \boldsymbol{u}|_{\partial\Omega} = \boldsymbol{u}_{0} \}, \]

或者等价地，使用应变表示的最小势能原理： \[ \begin{equation} \check\Pi = \min_{\boldsymbol{\varepsilon}\in \mathcal{E}} \Pi(\boldsymbol{\varepsilon}) = \min_{\boldsymbol{\varepsilon}\in \mathcal{E}}\frac{1}{2}\int_{\Omega} \boldsymbol{\varepsilon}:\mathbb{L}:\boldsymbol{\varepsilon}\ \mathrm{d}\Omega, \quad \mathcal{S} \triangleq \{ \boldsymbol{\varepsilon} \mid\boldsymbol{\varepsilon}=\nabla_{s}\boldsymbol{u},\ \boldsymbol{u} \in H^{1},\ \boldsymbol{u}|_{\partial\Omega} = \boldsymbol{u}_{0} \}, \label{eq:min_potl} \end{equation} \] 现在，引入参考模量 \(\mathbb{L}_{0}\)，将泛函 \(\eqref{eq:min_potl}\) 分解成如下两部分： \[ \Pi(\boldsymbol{\varepsilon}) = \underbrace{\frac{1}{2}\int_{\Omega} \boldsymbol{\varepsilon}:(\mathbb{L} - \mathbb{L}_{0}):\boldsymbol{\varepsilon} \ \mathrm{d}\Omega}_{\Delta \Pi} + \underbrace{\frac{1}{2}\int_{\Omega} \boldsymbol{\varepsilon}:\mathbb{L}_{0}:\boldsymbol{\varepsilon} \ \mathrm{d}\Omega}_{\Pi_{0}}. \] 假设 \(\Delta \mathbb{L} \triangleq \mathbb{L}-\mathbb{L}_{0}\) 正定，应用 Young-Fechel 变换，将 \(\Delta \Pi\) 表示为关于另一个二阶张量场 \(\boldsymbol{p}\) 的二次型变分问题： \[ \begin{equation} \Delta \Pi = \max_{\boldsymbol{p} \in \mathcal{P}} \left\{ \int_{\Omega} \left( \boldsymbol{p}:\boldsymbol{\varepsilon} - \frac{1}{2} \boldsymbol{p}:\Delta\mathbb{L}^{-1}:\boldsymbol{p} \right) \mathrm{d}\Omega \right\}, \quad \mathcal{P}= \{ \boldsymbol{p} \mid \boldsymbol{p} \in L^{2} \}. \label{eq:yf} \end{equation} \] 若 \(\Delta \mathbb{L} \triangleq \mathbb{L}-\mathbb{L}_{0}\) 负定，那么上式中求最大值将替换为求最小值。将式 \(\eqref{eq:yf}\) 代入泛函 \(\eqref{eq:min_potl}\) 中，并交换求最大最小的顺序，就得到关于两场泛函 \(\mathcal{H}(\boldsymbol{p},\boldsymbol{\varepsilon})\) 的变分问题： \[ \begin{equation} \check{\Pi} = \begin{cases} \max\limits_{\boldsymbol{p}\in\mathcal{P}} \min\limits_{\boldsymbol{\varepsilon}\in\mathcal{E}} \mathcal{H}(\boldsymbol{p},\boldsymbol{\varepsilon}), \quad \Delta\mathbb{L}\text{ is positive definite}, \\ \min\limits_{\boldsymbol{p}\in\mathcal{P}} \min\limits_{\boldsymbol{\varepsilon}\in\mathcal{E}} \mathcal{H}(\boldsymbol{p},\boldsymbol{\varepsilon}), \quad \Delta\mathbb{L}\text{ is negative definite}, \end{cases} \label{eq:hs_vari} \end{equation} \] 泛函 \(\mathcal{H}(\boldsymbol{p},\boldsymbol{\varepsilon})\) 的形式为： \[ \begin{equation} \mathcal{H}(\boldsymbol{p},\boldsymbol{\varepsilon}) \triangleq \int_{\Omega} \left( \boldsymbol{p}:\boldsymbol{\varepsilon} - \frac{1}{2} \boldsymbol{p}:\Delta\mathbb{L}^{-1}:\boldsymbol{p} + \frac{1}{2} \boldsymbol{\varepsilon}:\mathbb{L}_{0}:\boldsymbol{\varepsilon} \right) \mathrm{d}\Omega. \label{eq:hsfctl_ori} \end{equation} \] 接下来将通过分解位移场进一步简化泛函 \(\mathcal{H}\) 的形式。位移场 \(\boldsymbol{u}\) 可以分解为如下两部分：\(\boldsymbol{u} = \boldsymbol{u}_{0} + \tilde{\boldsymbol{u}}\)，其中 \(\boldsymbol{u}_{0}\) 是使得均质材料的势能 \(\Pi_{0}\) 取最小值的位移场，满足给定的位移边界条件，作为泛函中的已知量： \[ \check\Pi_{0}\triangleq\Pi_{0}(\boldsymbol{u}_{0}) = \min_{\boldsymbol{u}\in\mathcal{U}}\frac{1}{2}\int_{\Omega} \nabla_{s}\boldsymbol{u}:\mathbb{L}_{0}:\nabla_{s}\boldsymbol{u} \ \mathrm{d}\Omega; \] 而位移扰动场 \(\tilde{\boldsymbol{u}}\) 是新的依赖变量，在可行域空间 \(\mathcal{U}_{0}\triangleq \{ \boldsymbol{u} \mid \boldsymbol{u} \in H^{1},\ \boldsymbol{u}|_{\partial\Omega} = \boldsymbol{0} \}\) 中取值。相应的，应变场也可以分解为 \[ \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}_{0} + \tilde{\boldsymbol{\varepsilon}}, \quad \boldsymbol{\varepsilon}_{0} = \nabla_{s} \boldsymbol{u}_{0}, \quad \tilde{\boldsymbol{\varepsilon}} = \nabla_{s} \tilde{\boldsymbol{u}} \in \mathcal{S}_{0}. \] 按照如上分解，泛函 \(\eqref{eq:hsfctl_ori}\) 就可以简化为如下形式： \[ \begin{equation} \mathcal{H}(\boldsymbol{p},\tilde{\boldsymbol{\varepsilon}}) \triangleq \int_{\Omega} \left( \frac{1}{2}\tilde{\boldsymbol{\varepsilon}}:\mathbb{L}_{0}:\tilde{\boldsymbol{\varepsilon}} + \boldsymbol{p}:\tilde{\boldsymbol{\varepsilon}} + \boldsymbol{p}:\boldsymbol{\varepsilon}_{0} - \frac{1}{2} \boldsymbol{p}:\Delta\mathbb{L}^{-1}:\boldsymbol{p} \right) \mathrm{d}\Omega + \check\Pi_{0}, \label{eq:hsfctl} \end{equation} \] 而变分原理 \(\eqref{eq:hs_vari}\) 中应变场的可行域应改成 \(\mathcal{S}_{0}\)。

Hashin-Shtrikman 变分原理的不同表述

本文旨在阐述，为什么 Hashin-Shtrikman 变分原理能够得到相对 Voigt-Reuss 上下界更紧致的估计，并作为一种常见的变分技巧被广泛应用？一个简短的答案是：H-S 提供了一种技巧，放松了试验函数空间对连续性的约束条件。这就可以显式构造更复杂的试验函数，所以得到更准确的估计。

Voigt-Reuss 上下界

如果指定线性位移边界条件 \(\boldsymbol{u}|_{\partial\Omega} = \boldsymbol{\varepsilon}_{0}\cdot\boldsymbol{x}\)，那么最小势能原理要求真实位移解 \(\check{\boldsymbol{u}}\) 是如下变分问题的解： \[ \begin{equation} \check{\boldsymbol{u}} =\mathop{\mathrm{argmin}}_{\boldsymbol{u}\in \mathcal{U}}\ \Pi_{\boldsymbol{u}}(\boldsymbol{u}) = \frac{1}{2}\int_{\Omega} \nabla_{s}\boldsymbol{u}:\mathbb{L}:\nabla_{s}\boldsymbol{u} \ \mathrm{d}\Omega, \quad \mathcal{U} \triangleq \{ \boldsymbol{u} \mid \boldsymbol{u} \in H^{1},\ \boldsymbol{u}|_{\partial\Omega} = \boldsymbol{\varepsilon}_{0}\cdot\boldsymbol{x} \}, \label{eq:min_potential} \end{equation} \]

以及可以通过上述泛函的最小值 \(\check{\Pi}_{\boldsymbol{u}}\) 定义等效模量 \(\mathbb{L}_{\boldsymbol{u}}^{c}\)： \[ \frac{1}{2}|\Omega| \boldsymbol{\varepsilon}_{0} : \mathbb{L}_{\boldsymbol{u}}^{c} : \boldsymbol{\varepsilon}_{0} \triangleq \check{\Pi}_{\boldsymbol{u}}. \] 上下界的估计需要显示给出位移场的表达式，并且满足式 \(\eqref{eq:min_potential}\) 中的约束条件。最简单的构造方式是将位移场在边界处的表达式直接延拓到区域 \(\Omega\) 内，\(\hat{\boldsymbol{u}} = \boldsymbol{\varepsilon}_{0} \cdot\boldsymbol{x}\)，代入到变分问题中得到对等效模量上界的 Voigt 估计 \[ \frac{1}{2}|\Omega| \boldsymbol{\varepsilon}_{0} : \mathbb{L}_{\boldsymbol{u}}^{c} : \boldsymbol{\varepsilon}_{0} \leq \frac{1}{2}|\Omega| \boldsymbol{\varepsilon}_{0} : \left\langle \mathbb{L} \right\rangle : \boldsymbol{\varepsilon}_{0} \Rightarrow \mathbb{L}_{\boldsymbol{u}}^{c} \leq \mathbb{L}. \] 类似的，如果指定均匀应力边界条件 \(\boldsymbol{\sigma}|_{\partial \Omega} = \boldsymbol{\sigma}_{0}\)，那么最小余能原理要求真实的应力场 \(\check{\boldsymbol{\sigma}}\) 是如下变分问题的解： \[ \begin{equation} \check{\boldsymbol{\sigma}} =\mathop{\mathrm{argmin}}_{\boldsymbol{\sigma}\in \mathcal{S}}\ \Pi_{\boldsymbol{\sigma}}(\boldsymbol{\sigma}) = \frac{1}{2}\int_{\Omega} \boldsymbol{\sigma}:\mathbb{M}:\boldsymbol{\sigma} \ \mathrm{d}\Omega, \quad \mathcal{S} \triangleq \{ \boldsymbol{\sigma} \mid \boldsymbol{\sigma} \in L^{2},\ \boldsymbol{\sigma}|_{\partial \Omega} = \boldsymbol{\sigma}_{0},\ \boldsymbol{\sigma} \cdot \nabla = \boldsymbol{0} \}, \label{eq:min_cmpl} \end{equation} \] 并通过泛函的最小值 \(\check{\Pi}_{\boldsymbol{\sigma}}\) 定义等效模量 \(\mathbb{L}_{\boldsymbol{\sigma}}^{c}\)： \[ \frac{1}{2}|\Omega| \boldsymbol{\sigma}_{0} : \left( \mathbb{L}_{\boldsymbol{\sigma}}^{c} \right)^{-1} : \boldsymbol{\sigma}_{0} \triangleq \check{\Pi}_{\boldsymbol{\sigma}}. \] 同样，使用边界到区域内的自然延拓（也即内部应力场均匀分布），就可以得到对等效模量下界的 Reuss 估计： \[ \frac{1}{2}|\Omega| \boldsymbol{\sigma}_{0} : (\mathbb{L}_{\boldsymbol{\sigma}}^{c})^{-1} : \boldsymbol{\sigma}_{0} \leq \frac{1}{2}|\Omega| \boldsymbol{\sigma}_{0} : \langle \mathbb{L}^{-1} \rangle : \boldsymbol{\sigma}_{0} \Rightarrow \langle \mathbb{L}^{-1} \rangle^{-1} \leq \mathbb{L}_{\boldsymbol{\sigma}}^{c}. \] 注意，以上分别是在不同的边界条件下得到的等效模量 \(\mathbb{L}_{\boldsymbol{u}}^{c}\) 和 \(\mathbb{L}_{\boldsymbol{\sigma}}^{c}\)，一般并不相同，但如果应用 Hill 给出的代表性体积元的概念——等效模量与边界条件的类型无关——重定义与边界条件无关的等效模量 \(\mathbb{L}^{c}\triangleq\mathbb{L}_{\boldsymbol{u}}^{c}=\mathbb{L}_{\boldsymbol{\sigma}}^{c}\)，并有如下估计： \[ \begin{equation} \langle \mathbb{L}^{-1} \rangle^{-1} \leq \mathbb{L}^{c} \leq \langle \mathbb{L} \rangle. \label{eq:vrb} \end{equation} \] 如果考虑弹性模量是分片常值分布的，那么上式可以进一步化简为 \[ \begin{equation} \left( \sum_{\alpha=1}^{N} c^{(\alpha)} \mathbb{M}^{(\alpha)} \right)^{-1} \leq \mathbb{L}^{c} \leq \sum_{\alpha=1}^{N} c^{(\alpha)} \mathbb{L}^{(\alpha)} \label{eq:vrb_d} \end{equation} \]

应力场的连续性条件

考虑非均质的 RVE，包含任意多个夹杂相。RVE 体力项 \(\boldsymbol{f}=\boldsymbol{0}\)，边界条件为周期性边界条件。首先构造两组容许空间中的应力和应变场，分别为 \((\boldsymbol{\sigma}^{(1)}, \boldsymbol{\varepsilon}^{(1)})\) 和 \((\boldsymbol{\sigma}^{(2)}, \boldsymbol{\varepsilon}^{(2)})\)。第一组考虑温度均匀变化，由此在每一相上产生热应变： \[ \langle \boldsymbol{\sigma}^{(1)} \rangle=\boldsymbol{0}, \quad \langle \boldsymbol{\varepsilon}^{(1)} \rangle = \boldsymbol{m}^{c} \Delta\theta, \quad \boldsymbol{\varepsilon}^{(1)} = \mathbb{M}:\boldsymbol{\sigma}^{(1)} + \boldsymbol{m} \Delta \theta. \] 式中，\(\boldsymbol{m}\) 是材料的热膨胀二阶张量，在空间中可以非均匀分布；\(\boldsymbol{m}^{c}\) 是均匀化后得到的宏观热膨胀张量。第二组考虑在 RVE 上施加宏观应力 \(\boldsymbol{\sigma}_{0}\)，无温度变化，由此得到的单胞内应变与应力场的响应： \[ \langle \boldsymbol{\sigma}^{(2)} \rangle=\boldsymbol{\sigma}_{0}, \quad \langle \boldsymbol{\varepsilon} \rangle^{(2)} = \mathbb{M}^{c} : \boldsymbol{\sigma}_{0}, \quad \boldsymbol{\sigma}^{(2)} = \mathbb{F}:\boldsymbol{\sigma}_{0}, \quad \boldsymbol{\varepsilon}^{(2)} = \mathbb{M}:\mathbb{F}:\boldsymbol{\sigma}_{0}. \] 式中，四阶张量 \(\mathbb{F}\) 是应力影响函数，通过应力影响函数和宏观应力，可以还原细观尺度下的应力场。

在构造完满足平衡方程的应力场和满足相容条件的应变场之后，接下来将应用虚功原理，得到宏观热膨胀张量 \(\boldsymbol{m}^{c}\) 的表达式。注意到应变场 \(\boldsymbol{\varepsilon}^{(1)}\) 和 \(\boldsymbol{\varepsilon}^{(2)}\) 体积平均后不等于零，所以不能用做虚位移表示的虚功原理中应变场的扰动，因此使用由虚应力表示的虚功原理： \[ \begin{equation} \int_{\Omega} \delta\boldsymbol{\sigma} : \boldsymbol{\varepsilon} \,\mathrm{d}V - |\Omega| \langle \delta\boldsymbol{\sigma} \rangle: \langle \boldsymbol{\varepsilon} \rangle = 0, \label{eq:vs} \end{equation} \] 其中，虚应力要求满足平衡方程 \(\delta\boldsymbol{\sigma}\cdot\nabla=\boldsymbol{0}\)，以及周期性边界条件；应变场只需要满足相容条件即可。首先将 \(\boldsymbol{\sigma}^{(1)}\) 作为虚应力，\(\boldsymbol{\varepsilon}^{(2)}\) 作为试验函数，代入虚功原理 \(\eqref{eq:vs}\) 中得到： \[ \int_{\Omega} \boldsymbol{\sigma}^{(1)} : \mathbb{M} : \boldsymbol{\sigma}^{(2)} \,\mathrm{d}V = 0. \] 然后再将 \(\boldsymbol{\sigma}^{(2)}\) 和 \(\boldsymbol{\varepsilon}^{(1)}\) 代入得到 \[ |\Omega|\boldsymbol{\sigma}_{0} : \boldsymbol{m}^{c}\Delta\theta = \int_{\Omega} (\mathbb{F}:\boldsymbol{\sigma}_{0}) : \mathbb{M}:\boldsymbol{\sigma}^{(1)} \,\mathrm{d}V + \Delta\theta\int_{\Omega} (\mathbb{F}:\boldsymbol{\sigma}_{0}) : \boldsymbol{m} \,\mathrm{d}V. \] 因为 \(\boldsymbol{\sigma}^{(1)}\) 的体积平均等于零，所以上式等号右侧第一项积分等于零（这一结论同样来源于 \(\eqref{eq:vs}\)）。再根据 \(\boldsymbol{\sigma}_{0}\) 的任意性，就得到

\[ \boldsymbol{m}^c = \frac{1}{|\Omega|} \int_{\Omega} \boldsymbol{m}: \mathbb{F}\,\mathrm{d}V. \] 如果令 \(\boldsymbol{\mu} = \boldsymbol{m} \Delta \theta\)，\(\boldsymbol{\mu}\) 是本征应变，那么上式就可以写作 \[ \begin{equation} \boxed{\boldsymbol{\mu}^c = \frac{1}{|\Omega|} \int_{\Omega} \boldsymbol{\mu}: \mathbb{F}\,\mathrm{d}V.} \label{eq:levin} \end{equation} \] 进一步的，如果本征应变场分片常值分布，那么上式又可以写成求和的形式： \[ \begin{equation} \boxed{\boldsymbol{\mu}^c = \sum_{\alpha=1}^{N} c^{(\alpha)} \boldsymbol{\boldsymbol{\mu}}^{(\alpha)} : \mathbb{F}^{(\alpha)}.} \label{eq:levin_sum} \end{equation} \] 式中，\(c^{(\alpha)}\) 是体积分数，\(\boldsymbol{\mu}^{(\alpha)}\) 是区域 \(\Omega^{(\alpha)}\) 内常值分布的本征应变，\(\mathbb{F}^{(\alpha)}\) 是在区域 \(\Omega^{(\alpha)}\) 内对应力影响函数的体积平均。

Remark

• 宏观尺度的应力和应变等于细观尺度中应力和应变场的体积平均，因此符号 \(\boldsymbol{\sigma}^{c}\) 和 \(\langle \boldsymbol{\sigma} \rangle\)，\(\boldsymbol{\varepsilon}^{c}\) 和 \(\langle \boldsymbol{\varepsilon} \rangle\) 是可以不加区别地混用的；然而，宏观本征应变并不是本征应变场的体积平均，而是需要应力影响函数进行加权，如公式 \(\eqref{eq:levin}\) 和 \(\eqref{eq:levin_sum}\) 所示。
• 公式的推导过程只用到虚功原理，因此对于任意非弹性变形，在用本征应变进行表示之后，Levin 总是成立。

这篇文章首先给出 Gâteaux 导数 与 Fréchet 导数的定义和符号，然后提供三个简单算例：（1）一维空间中的导数；（2）二维空间的方向导数，并给出 Fréchet 导数的图像解释；（3）无穷维空间中泛函的变分。

Gâteaux 导数 与 Fréchet 导数的定义

Gâteaux 导数的定义为 \[ D_{g} F( y;\eta ) = \lim_{\epsilon\to0} \frac{F( y+\epsilon \eta ) - F( y)}{\epsilon} = \left.\frac{\mathrm{d} }{\mathrm{d} \epsilon} F( y+\epsilon \eta )\right|_{\epsilon=0} \] 由定义式可看到，如果将给定点处函数 \(F\) 的 Gâteaux 导数看作是关于 \(\eta\) 的新函数，那么 \(D_{g}F\) 与 \(F\) 属于同一个函数空间。当涉及到定义在无穷维函数空间中的泛函时，需要严格界定泛函定义域所属的函数空间，以及测试函数 \(\eta\) 所在的函数空间。

此外，根据定义式，Gâteaux 导数是关于函数 \(\eta\) 的一阶齐次函数，也即 \[ D_{g}F(y;\alpha \eta) = \alpha D_{g}F(y;\eta), \quad \forall \alpha \in \R \] 继续推广 Gâteaux 导数的定义到高阶情景，有 \[ D_{g}^{n} F( y;\eta ) = \left.\frac{\mathrm{d}^n }{\mathrm{d} \epsilon^n} F( y+\epsilon \eta )\right|_{\epsilon=0} \]

Fréchet 导数是在点 \(y\) 处寻找一个有界线性映射 \(D_{F}: X\mapsto Y\)，其中 \(X\) 是被求导函数定义域所属的空间，\(Y\) 是被求导函数值域所属的空间。若关于点 \(y\)，对任意的 \(h\in X\)，存在有界线性映射 \(D_{F}(y)\)，使得 \[ \lim_{\|h\|\to 0} \frac{ \| F(y+h) - F(y) -D_{F}(y)\cdot h \| }{\|h\|} = 0 \] 那么函数在该点处 Fréchet 可微。

这篇文档整理了力学中常用的能量原理在三维空间小变形问题背景下的一般性表述，包含（1）容许空间的概念（2）虚功（虚位移/虚应力）原理，以及限制在弹性问题背景下的（3）最小势能和最小余能原理（4）Clapeyron 定理，和（5）Betti 互易定理。为使得这篇文档不会变得过分冗长，作者默认读者熟悉基本的变分概念和变分的计算方法，只给出能量原理的推导过程，没有给出具体的应用实例。

参考文献

• Reddy, J N. Energy Principles and Variational Methods in Applied Mechanics.
• Dvorak, George. Micromechanics of Composite Materials.

容许空间

首先，这些原理的问题背景是区域 \(\Omega\) 上的边值问题（boundary value problem, BVP）。其边界条件为力边界和位移边界的组合，\(\partial \Omega = \partial\Omega_{t} \cup \partial\Omega_{u}\)，位移边界条件为 \(\boldsymbol{u}=\boldsymbol{u}_0\)，力边界条件为 \(\boldsymbol{t}=\boldsymbol{\sigma}\cdot\boldsymbol{n}\)。区域内满足平衡方程和边界条件： \[ \boldsymbol{\sigma} \cdot \nabla + \boldsymbol{f} = \boldsymbol{0},\quad \boldsymbol{\sigma}\cdot\boldsymbol{n}=\boldsymbol{t}_{0} \text{ at } \partial_{t}\Omega, \quad \boldsymbol{u}=\boldsymbol{u}_{0} \text{ at } \partial_{u}\Omega. \] 应变场满足协调方程： \[ \boldsymbol{\varepsilon} = \frac{1}{2} ( \boldsymbol{u}\nabla + \nabla\boldsymbol{u} )\triangleq\nabla_{s}\boldsymbol{u}. \] 应力和应变之间由本构方程关联。在上述问题背景下，一些常见于原理的概念列举如下：

这篇文章旨在解释为什么在计算变分时，在形式上与计算微分是如此相似。这需要在符号以及符号遵循的规则之外，理解微分或变分究竟是什么。如果将一阶微分或变分解释成各自对偶空间中的线性泛函，那么就可以通过线性泛函的运算规律，说明为什么微分和变分运算在形式上是一致的。

一阶变分的重定义

对如下泛函 \[ \begin{equation} \Pi(\boldsymbol{u}) = \int_{\Omega} F(\boldsymbol{x},\boldsymbol{u},\nabla\boldsymbol{u}) \,\mathrm{d}V, \quad \boldsymbol{u}\in \mathcal{V}, \quad \mathcal{V} = \{ \boldsymbol{v} \mid \boldsymbol{v}\in H^1, \boldsymbol{v}|_{\partial\Omega_{u}} = \boldsymbol{v}_{0} \}, \label{eq:functional_def} \end{equation} \] 其关于 \(\boldsymbol{\eta}\) 的 Gâteaux 导数等于 \[ \begin{equation} D_{g}\Pi(\boldsymbol{u};\boldsymbol{\eta}) = \int_{\Omega} \frac{\partial F}{\partial \boldsymbol{u}} \cdot \boldsymbol{\eta} + \frac{\partial F}{\partial \nabla\boldsymbol{u}}:\nabla\boldsymbol{\eta} \,\mathrm{d}V. \label{eq:vari_1} \end{equation} \]

容易验证，式 \(\eqref{eq:vari_1}\) 在点 \(\boldsymbol{u}\) 处定义了测试函数空间 \(\mathcal{V}_{0}=\{ \boldsymbol{v} \mid \boldsymbol{v}\in H^1, \boldsymbol{v}|_{\partial\Omega_{u}} = \boldsymbol{0} \}\) 上的线性泛函，线性泛函组成的空间称为对偶空间，记作 \(\mathcal{V}_{0}^{*}\)。另一方面，根据变分的定义，该线性泛函其实就是泛函 \(\Pi\) 在点 \(\boldsymbol{u}\) 处的一阶变分： \[ \begin{equation} \delta\Pi(\boldsymbol{u}) = \int_{\Omega} \frac{\partial F}{\partial \boldsymbol{u}} \cdot \delta\boldsymbol{u} + \frac{\partial F}{\partial \nabla\boldsymbol{u}}:\nabla\delta\boldsymbol{u} \,\mathrm{d}V = D_{g}\Pi(\boldsymbol{u};\delta\boldsymbol{u}). \label{eq:vari_def} \end{equation} \]

因此，泛函在点 \(\boldsymbol{u}^{*}\) 处取极值的必要条件——其一阶变分等于零：\(\delta\Pi =0\)——另一种表述是：泛函在点 \(\boldsymbol{u}^{*}\) 处的 Gâteaux 导数是对偶空间 \(\mathcal{V}_{0}^{*}\) 中的零元素。

有限元离散后, 要求节点处内力与外力平衡, \(\boldsymbol{\mathsf{f}}^{\mathrm{int}}=\boldsymbol{\mathsf{f}}^{\mathrm{ext}}\). 其中, 内力可由刚度矩阵 \(\boldsymbol{\mathsf{K}}\) 和节点自由度列阵 \(\boldsymbol{\mathsf{d}}\) 给出, 最终得到离散后的平衡方程: \[ \begin{equation} \boldsymbol{\mathsf{K}} \boldsymbol{\mathsf{d}} = \boldsymbol{\mathsf{f}}^{\mathrm{ext}}. \end{equation} \] 除平衡方程之外, 节点之间一般还需满足若干运动学约束: \[ F_{k}( \boldsymbol{d}_{a}, \ldots ) = \boldsymbol{0}, \quad k=1,\ldots \]

这篇文档记录作者常用且反复查找的与 Abaqus 相关的命令和操作。

Pull-back 和 push-forword 操作

变形梯度不仅将参考构型中的微元 \(\mathrm{d}X\) 映射为当前构型中的微元 \(\mathrm{d}x\), 同时作为两点张量, 还将物质形式的物理量映射到空间形式的物理量. 因此, 对于空间向量 \(v_{i}(x,t)\), 可以从数学的意义上, 将该变量 pull back 到物质坐标表示的物质形式的向量 \(V_I(X,t)\): \[ V_I(X,t) \triangleq F_{Ii}^{-1} v_{i}(\varphi(X,t), t) \] 类似的, 也可以将用物质形式的向量 \(G_{I}\), push forward 到空间形式的向量 \(g_i\): \[ g_i(x,t) \triangleq F_{iI} G_{I}(\varphi^{-1}(x,t), t) \] 如果使用变形梯度作为对向量前推后拉的操作运算, 那么除了使用 \(F_{iI}\), \(F_{Ii}^{-1}\) 之外, 还可以使用变形梯度的转置 \(F_{iI}^{-\top}\), \(F_{Ii}^{\top}\). 对每一个指标, 有两种前推 \(\{ F_{iI}, F_{iI}^{-\top} \}\) 后拉 \(\{ F_{Ii}^{-1}, F_{Ii}^{T} \}\) 的方式. 因此, 如果推广到对二阶张量的前推后拉操作, 每种运算有四种方式.

Cauchy 应力的定义

在 Spencer 的连续介质力学书中, Cauchy 应力定义为 \[ t^{(i)} = \sigma_{ij} e_{j}, \quad \mathrm{or} \quad \sigma_{ij} \triangleq t^{(i)} \cdot e_{j} \] 这样得到的平衡方程为 \[ \nabla \cdot \sigma + \rho b = 0, \quad \mathrm{or} \quad \sigma_{ji,j} + \rho b_{i} = 0 \tag{1} \] 而我熟悉的平衡方程为 \[ \sigma_{ij,j} + \rho b_i = 0 \]

空间坐标下的守恒方程可通过输运定理获得, 而对于物质坐标下的守恒方程, 通过积分参数变换的方式获得. 因此在本文中, 首先陈述微元在不同参数表示下的变换关系, 然后再根据空间坐标下的守恒方程, 通过参数变换, 得到物质坐标下的守恒方程.

对于 \(n\) 阶行列式 \(|A|\), 去掉 \(i\) 行和 \(j\) 列得到的 \(n-1\) 阶行列式记作余因式 \(c_{ij}\), 行列式 \(|A|\) 的代数余子式定义为 \[ m_{ij} \triangleq (-1)^{i+j} c_{ij} \] 行列式可以按照行或者列展开为低阶行列式的线性组合: \[ |A| = \sum_{i=1}^{n} a_{ij} m_{ij} = \sum_{j=1}^{n} a_{ij} m_{ij} \]

质量守恒

考虑空间坐标下的控制体 \(V\), 在 \(\delta t\) 时间内, 流出控制体 \(V\) 的质量为 \[ \mathrm{d} m = \delta_t\int_{\partial V} \rho v\cdot n ~\mathrm{d}S + \mathcal{O}(\delta t^2) \] 控制体在 \(t+\delta t\) 时刻的质量为 \[ m_V(t+\delta t) = \int_{V} \rho(x,t+\delta t) ~\mathrm{d}x \] 质量守恒的表述为: 在控制体的质量改变, 只能通过质量进出控制体引起. 因此有 \[ m_{V}(t + \delta t) = m_{V}(t) - \mathrm{d} m \Rightarrow m_{V}(t + \delta t) - m_{V}(t) + \mathrm{d} m = 0 \]

小变形问题

在有限变形的框架下, 使用物质坐标系和空间坐标系描述应变的度量分别为 Lagrange 应变和 Euler 应变: \[ E \triangleq \frac{1}{2}(F^{\top} F - I), \quad e \triangleq \frac{1}{2}(I - F^{-\top} F^{-1}) \tag{1} \] 式中的 \(F\) 和 \(F^{-1}\) 分别为 \[ F_{iI} = \frac{\partial x_i}{\partial X_I}, \quad F_{Ii}^{-1} = \frac{\partial X_I}{\partial x_i} \]

向量空间的定义

向量空间的定义在此不再赘述, 不过需要补充一些 remarks:

[!NOTE]

1. 向量空间只定义了运算的概念, 并没有引入距离, 角度等拓扑概念.
2. 向量空间中的“向量”并不局限于后文提到的一阶张量, 而是更抽象的概念, 比如表示全体连续函数的向量空间 \(C^{0}\), 连续函数 \(f\) 是向量空间 \(C^0\) 中的一个元素.
3. 向量空间总是需要伴随一个数域 \(\mathbb{F}\), 可以是实数域 \(\mathbb{R}\), 也可以是复数域 \(\mathbb{C}\). 对于有限维的向量空间, 在给定一组基底之后, 就能够定义一个从 \(V\) 到 \(\mathbb{F}^n\) (\(n\) 为向量空间的维度) 的一一映射, 因此 \(V\) 同构于 \(\mathbb{F}^n\)

首先, 先定义梯度, 散度和旋度的概念, 并特别强调对于高阶张量的定义, 这将规范使用的 notation. 然后, 以二维空间为例, 引出散度定理, Stokes 公式和 Gauss 公式. 最后, 再给出一般的张量场中这三种公式的形式.

梯度, 散度和旋度

引入记号 \(\nabla\), \(\nabla \cdot\), \(\nabla \times\), 分别表示梯度, 散度和旋度. 首先介绍这三种运算最简单的情况, 然后再推广到一般的张量空间当中.

对标量场的梯度运算, 得到的是一个向量场: \[ \nabla \phi = \frac{\partial \phi\ }{\partial x_i} e_i \] 对向量场的散度运算, 得到的是一个标量场: \[ \nabla \cdot v = \frac{\partial v_i}{\partial x_i} \] 对向量场的旋度运算, 得到的仍是一个向量场: \[ \nabla \times v = \begin{vmatrix} i & j & k \\ \dfrac{\partial \ }{\partial x_1} & \dfrac{\partial \ }{\partial x_1} & \dfrac{\partial \ }{\partial x_1} \\ v_1 & v_2 & v_3 \end{vmatrix} \]