有限元中的能量原理
初边值问题的陈述
物质描述
使用物质描述的控制方程定义在 \(B_0 \times (0,t_0)\) 区域上, 方程的未知量为 \(\varphi(X,t)\), 满足守恒定律的控制方程为 \[ (\varphi \nabla_0) \cdot S \cdot \nabla_0 + \rho_0 b = \rho_0 \frac{ \partial^2 \varphi }{ \partial t^2 } \tag{1.1.a} \] 式中, \(S\) 为 PK2 应力. 对于是初边值问题, 定解需要给定初始条件和边界条件. 其中初始条件为 \[ \varphi(X,0) = X, \quad \frac{ \partial \varphi }{ \partial t}(X,0) = v_0(X) \tag{1.1.b} \] 方程的边界条件有两类, 一类是位移边界条件, 在参考构型中的边界处指定 \(t>0\) 时刻的位移: \[ \varphi(X,t) = \overline{x}(X,t), \quad X \in \partial B_{0u},~ t > 0 \tag{1.1.c} \] 另一类是力边界条件, 在参考构型的边界处指定 \(t>0\) 时刻的面力: \[ F \cdot S \cdot N(X) = \overline{T}(X,t) \quad X \in \partial B_{0t},~ t > 0 \tag{1.1.d} \] 这两类边界条件不能重叠, 也即不能在边界点同一个分量处指定两种边界条件, 否则就是 ill-posed: \[ \partial B_{0u} \cup \partial B_{0t} = \partial B_{0}, \quad \partial B_{0u} \cap \partial B_{0t} = \empty \] 除上述初边值条件之外, 还需要提供本构关系, 将 PK2 应力与变形和变形历史相关联, 以及 \(t=0\) 时刻的密度 \(\rho_0 = \rho_0(X)\) 和 \(t>0\) 时刻的体积力 \(b = b(X,t)\). 至此, 待求解的力学初边值问题 (1) 才算完整地陈述.
[!NOTE]
使用位移 \(u\) 作为未知量表述的初边值问题
方程中的未知量也可以是位移 \(u\), 与 \(\varphi\) 的关系为 \[ u(X,t) = \varphi(X,t) - X \tag{1.2} \] 相应的初始条件为 \[ u(X,0) = 0, \quad \frac{ \partial u }{ \partial t}(X,0) = v_0(X) \] 以及位移边界条件为 \[ u(X,t) = \overline{u}(X,t), \quad X \in \partial B_{0u},~ t > 0 \] 上式与式 (1.1.c) 的关系为 \[ \overline{u}(X,t) = \overline{x}(X,t) - X \]
空间描述
势能原理
在这一节中, 我们将平衡方程与某一能量泛函取驻值的条件等价, 由此得到势能驻值原理. 泛函的驻值点可能是极大值或极小值或鞍点, 我们将指出系统的稳定平衡状态对应能量泛函的极小值, 这也是最小势能原理.
平衡方程的势能泛函
若不考虑方程 (1.1.a) 中的加速度项, 就得到材料在平衡状态下满足的控制方程为 \[ (\varphi \nabla_0) \cdot S \cdot \nabla_0 + \rho_0 b = 0 \tag{2.1.a} \] 以及边界条件 \[ \varphi(X) = \overline{x}(X), \quad X \in \partial B_{0u} \\ F \cdot S \cdot N(X) = \overline{T}(X) \quad X \in \partial B_{0t} \tag{2.1.b} \] 接下来, 我们将构造一个与上述边值问题等价的变分问题. 首先构造如下能量泛函 \(\Pi:L \mapsto \mathbb{R}\), 等于材料中的应变能与外力势能之和: \[ \Pi(\varphi) = \underbrace{\int_{B_0} W(F) ~\mathrm{d} X }_{应变能} \underbrace{ - \int_{B_0} \rho_0 b \cdot \varphi ~\mathrm{d} X - \int_{\partial B_{0t}} \overline{T} \cdot \varphi ~\mathrm{d} \Sigma(X) }_{外力势能} \tag{2.2.a} \] 式中, \(W\) 为应变能密度函数 (质量密度), 变形梯度 \(F = \varphi \nabla_0\). 泛函 \(\Pi\) 定义在函数空间上, 为了建立等价关系, 除了确定泛函的形式之外, 还需要确定函数变量 \(\varphi\) 所在的空间 \(L\). 我们约束变形映射 \(\varphi\) 满足与边值问题的位移边界条件 (2.1.b) 一致: \[ \varphi \in L = \{ f \mid f \in C^0(B_0), \quad f = \overline{x} ~~\text{on}~~ \partial B_{0u} \} \tag{2.2.b} \]
泛函的变分
在函数空间 \(L\) 中点 \(\varphi_0\) 附近, 定义一条经过点 \(\varphi_0\) 的参数曲线: \[ \varphi(\varepsilon) = \varphi_0 + \varepsilon \eta, \quad \varepsilon \in \mathbb{R} \tag{2.3} \] 由此得到关于参数 \(\varepsilon\) 的函数 \(\phi(\varepsilon)\): \[ \phi(\varepsilon) \triangleq \Pi(\varphi_0 + \varepsilon \eta) \] 若泛函在点 \(\varphi_0\) 处取极大(小)值, 那么对应函数 \(\phi\) 在 \(\varepsilon = 0\) 处取极大(小)值, 因此 \[ \left. \frac{\mathrm{d} \phi}{\mathrm{d} \varepsilon} \right|_{\varepsilon=0} = 0 \tag{2.4} \] 是泛函取极值的一个必要条件. 如果定义运算 \(\delta_{\eta}: L \mapsto L\) \[ \boxed{ \delta_{\eta}[f](\varphi_0) \triangleq \left. \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} f(\varphi_0 + \varepsilon \eta) \right|_{\varepsilon=0} } \tag{2.5} \] 那么条件 (2.4) 就可以重新叙述为 \[ \delta_{\eta}[\Pi](\varphi_0) = 0, \quad \forall \eta \in \overline{L} \tag{2.6} \] 接下来确定 \(\eta\) 的取值空间. 根据式 (2.3), \(\varphi_0 \in L\), \(\varphi(\varepsilon) \in L\), 因此 \(\eta\) 的取值函数空间对应为位移边界处恒等于 0 的函数: \[ \eta \in \overline{L} = \{ f \mid f \in C^0(B_0), \quad f = 0 ~~\text{on}~~ \partial B_{0u} \} \tag{2.7} \] 为了在记号上与微分符号进行类比, 我们将 \(\varepsilon\eta\) 重映射为 \(\varepsilon\delta \varphi\), 表示关于独立变量 \(\varphi\) 的变分. 当独立变量的变分意义明确时, 省略运算符 \(\delta_{\eta}\) 的下标, 并将 \(\delta_{\eta}[f]\) 简记作 \(\delta f\). 根据式 (2.5), 运算符 \(\delta\) 的定义是基于单变量微积分的导数, 因此继承了关于求导数的线性性质, 所以可以将 \(\delta\) 与微分符号 \(~\mathrm{d}\) 进行类推.
[!TIP]
变分运算公式
函数的乘积: \[ \begin{aligned} \delta(fg) &= \left. \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} f(\varphi + \varepsilon \delta\varphi) g(\varphi + \varepsilon \delta\varphi) \right|_{\varepsilon=0} \\ &= \left( \left. f(\varphi + \varepsilon \delta\varphi) \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} g(\varphi + \varepsilon \delta\varphi) + g(\varphi + \varepsilon \delta\varphi) \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} f(\varphi + \varepsilon \delta\varphi) \right)\right|_{\varepsilon=0} \\ &= f ~\delta g + g ~\delta f \end{aligned} \] 微分运算: \[ \delta (\varphi \circ \nabla_0) = \left. \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} \left( (\varphi \circ \nabla_0) + \varepsilon (\delta \varphi \circ \nabla_0) \right) \right|_{\varepsilon=0} = \delta \varphi \circ \nabla_0 \] 式中, \(\circ\) 表示梯度, 散度 \(\cdot\), 叉乘 \(\times\) 中的一种.
积分运算: \[ \delta \int f(\varphi) ~\mathrm{d}X = \left. \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} \int f(\varphi + \varepsilon \delta \varphi) ~\mathrm{d}X \right|_{\varepsilon=0} = \int \left. \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} f(\varphi + \varepsilon \delta \varphi) \right|_{\varepsilon=0} ~\mathrm{d}X = \int \delta f ~\mathrm{d} X \] 复合函数: \[ \left. \delta f \circ g = \frac{\mathrm{d} ~}{\mathrm{d} \varepsilon} f \circ g (\varphi + \varepsilon \delta \varphi) \right|_{\varepsilon=0} = \frac{\mathrm{d} f}{\mathrm{d} g} \delta g \] 最后, 我们说明即使选取不同的独立变量 \(u\) 或者 \(\varphi\), 它们的变分总是相同的. 根据式 (1.2): \[ \delta u = \delta (\varphi -X) \] \(X\) 与 \(\varphi\) 无关, 因此有 \[ \delta u = \delta \varphi \] 因此在下文的推导中, 我们用更符合惯例的 \(\delta u\) 代替 \(\delta \varphi\).
泛函驻值条件蕴含平衡方程
根据式 (6), 泛函 \(\Pi\) 取驻值的必要条件为 \[ \delta \Pi = 0, \quad \forall \delta u \in \overline{L} \] 将式 (2.2.a) 代入上式中得到 \[ \delta \Pi = \delta \int_{B_0} W(F) ~\mathrm{d} X - \delta \int_{B_0} \rho_0 b \cdot \varphi ~\mathrm{d} X - \delta \int_{\partial B_{0t}} \overline{T} \cdot \varphi ~\mathrm{d} \Sigma(X) \equiv 0, \quad \forall \delta u \in \overline{L} \tag{2.8} \] 式中, \[ \delta \int_{B_0} W(F) ~\mathrm{d} X = \int_{B_0} \delta W(F) ~\mathrm{d} X = \int_{B_0} \frac{ \partial W }{ \partial F} : (\delta u \nabla_0) ~\mathrm{d} X \tag{2.9} \] 如果考虑材料为超弹材料, 那么 \[ P = \frac{ \partial W }{ \partial F} \] 式中, \(P\) 为 PK1 应力, \(P = F \cdot S\). 应用分部积分公式得到 \[ P : (\delta u \nabla_0) = (\delta u \cdot P) \cdot \nabla_0 - (P \cdot \nabla_0) \cdot \delta u \] 代入到 (2.9) 式右端, 并注意到变分 \(\delta u \equiv 0\) 在边界 \(\partial B_{0u}\): \[ \boxed{ \int_{B_0} P : (\delta u \nabla_0) ~\mathrm{d} X = \int_{\partial B_{0u}} \delta u \cdot P \cdot N ~\mathrm{d} \Sigma(X) - \int_{B_0} (P \cdot \nabla_0) \cdot \delta u ~\mathrm{d} X } \tag{2.10} \] 式 (2.10) 是虚功原理的特殊形式, 这里我们是从守恒系统的变分等于 0 出发, 得到式 (2.9). 之后我们将从平衡方程出发, 得到更一般形式的虚功原理.
将 (2.10) 代入到 (2.8) 中得到 \[ - \int_{B_0} (P \cdot \nabla_0 + \rho_0 b) \cdot \delta u ~\mathrm{d} X + \int_{\partial B_{0t}} (P \cdot N -\overline{T}) \cdot \delta u ~\mathrm{d} \Sigma(X) = 0, \quad \forall \delta u \in \overline{L} \] 再根据 \(\delta u\) 的任意性, 就得到平衡方程 (2.1.a) 和力边界条件 (2.1.b).
平衡态的稳定性, 最小势能原理
在上述推导中, 我们只预设了变分的驻值性, 并由此得到平衡方程和力边界条件. 进一步地, 如果我们能确定泛函在驻值点处如果取得极小值, 那么就可以确定系统的平衡态是稳定的. 这就是最小势能原理.
虚功原理
首先我们定义什么是 statically admissible forces and stresses 和 kinematically admissible displacements:
- statically admissible stresses: 满足平衡方程和力边界条件的应力场分布 \(\tilde P\);
- kinematically admissible displacements 满足位移边界条件的位移场分布.
在假想的位移场下, 可允许位移之间的差值被定义为虚位移场 \(\delta \tilde u\), 虚位移场在式 (2.7) 定义的函数空间当中. 当施加虚位移场 \(\delta \tilde u\) 时, 假设应力场和外力固定不变, 因此外力在虚位移场下所做的虚外力功为力边界处的面力与体积力做功之和: \[ \delta W_{E} \triangleq \int_{\partial B_{0t}} \delta \tilde u \cdot \overline{T} ~\mathrm{d}X + \int_{B_0} \rho_0 b \cdot \delta \tilde u ~\mathrm{d}\Sigma(X) \] 虚应力场 \(\tilde{P}\) 满足指定的力边界条件 \(\tilde{P} \cdot N = \overline{T}\), 因此有 \[ \int_{\partial B_{0t}} \delta \tilde u \cdot \overline{T} ~\mathrm{d}X = \int_{\partial B_{0t}} \delta \tilde u \cdot \tilde{P} \cdot N ~\mathrm{d}X \] 注意到 \(\delta \tilde u\) 在边界 \(\partial B_{0u}\) 处恒等于 0, 因此可以将上式中边界积分项拓展到封闭边界区域的积分. 再应用散度定理得到 \[ \oint_{\partial B_{0}} \delta \tilde u \cdot \tilde{P} \cdot N ~\mathrm{d}X = \int_{B_0} (\delta \tilde u \nabla_0) : \tilde{P} + \delta \tilde u \cdot (\tilde{P} \cdot \nabla_0) ~\mathrm{d}X \] 定义上式中第二项为内力虚功: \[ \delta W_{I} \triangleq \int_{B_0} (\delta \tilde u \nabla_0) : \tilde{P} \] 同时注意到 \(\tilde{P}\) 是满足平衡方程的虚应力场: \[ \tilde{P} \cdot \nabla_0 + \rho_0 b = 0 \] 最终我们得到 \[ \delta W_{E} = \delta W_{I} \] 因此, 对于满足平衡方程和力边界条件的虚应力场, 以及虚位移场, 外力对系统的虚位移所做的外力功, 等于虚应力场在虚位移导致的虚应变场所做的内力功. 在虚功原理的陈述中, 应力场和位移场是完全独立的, 没有考虑任何附加的本构关系.