虚功原理和最小势能原理推导

JCHEN161

已于 2024-11-12 00:30:55 修改

阅读量1.4k

点赞数 26

文章标签：学习

于 2024-11-11 09:11:04 首次发布

本文链接：https://blog.csdn.net/mrforth/article/details/143669472

版权

0. 写在前面

本文是在作者学习了部分弹性力学和变分的知识下做的一个笔记，其实对于标题所表明的推导我不能说已经完全掌握。这篇文章的起始点完全是由于我不满足于应用这两个原理，所以一直想知道如何证明和推导，然后发现还要学变分的内容，又懵懵懂懂地学了一段时间，直到现在也不能说学的多明白，只能用这篇笔记来记录下目前的成果。

1. 预备知识

首先我们定义如下概念：
面力：分布在物体表面的作用力，如摩擦力、接触力等

体力：分布在物体体积内部的力，如重力、磁力等

描述一个点的应力状态有：

$\sigma=\begin{pmatrix} \sigma_x & \sigma_y & \sigma_z & \tau_{yz}& \tau_{zx} & \tau_{xy} \end{pmatrix}^T$

算子矩阵：

$E(\nabla)=\begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 & 0 & \frac{\partial}{\partial z} &\frac{\partial}{\partial y} \\ 0 & \frac{\partial}{\partial y} & 0 & \frac{\partial}{\partial z} & 0 &\frac{\partial}{\partial x} \\ 0& 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{bmatrix}$

矩阵：

$E(n)=\begin{bmatrix} n_x & 0 & 0 & 0 & n_z & n_y\\ 0 & n_y & 0 & n_z & 0 & n_x\\ 0 & 0 & n_z & n_y & n_x & 0 \end{bmatrix}$

变形协调关系：

$\epsilon=E^T(\nabla)u$

$u=(u_x,u_y,u_z)^T$

平衡方程：

$E(\nabla)\sigma+f=0$

边界条件：

$E(n)\sigma = \bar p\\ u = \bar{u}$

2. 虚功原理

考虑一个连续体受到了面力和体力产生了变形，在该物体内部肯定产生了相应的应力和应变，那么考虑一个微分的体积 $d\Omega$ ，它的应变能为:

$\sigma^T\epsilon d\Omega=\sigma^T E^T(\nabla)ud\Omega$

该等式成立是根据变形协调关系，因此在整个体积上应变能：
$\iiint _{\Omega}\sigma^T E^T(\nabla)ud\Omega$

除了应变能之外再分别考虑面力和体力所做的功产生的外力势能。考虑一个微分的表面 $dS$ ，根据边界条件有：

$E(n)\sigma = \bar p$

所以面力所做的功为：

$(E(n)\sigma)^T udS$

在微分体积 $d\Omega$ 上，根据平衡方程我们知道

$f=-E(\nabla)\sigma$

所以体力所做功的大小为：

$(E(\nabla)\sigma)^T ud\Omega$

因此整体分别为：
$\iint_{\partial \Omega}(E(n)\sigma)^T udS$

$\iiint_{\Omega}(E(\nabla)\sigma)^T ud\Omega$

下面我们需要证明这个等式，即

$\iiint _{\Omega}\sigma^T E^T(\nabla)ud\Omega=\iint_{\partial \Omega}(E(n)\sigma)^T udS-\iiint_{\Omega}(E(\nabla)\sigma)^T ud\Omega$

从纯数学的角度出发，首先规定

$E(\nabla)=E_1\frac{\partial}{\partial x}+E_2\frac{\partial}{\partial y}+E_3\frac{\partial}{\partial z}$

$(E(\nabla)\sigma)^T u=(\frac{\partial \sigma^T}{\partial x}E_1^T+\frac{\partial \sigma^T}{\partial y}E_2^T+\frac{\partial \sigma^T}{\partial z}E_3^T)u$

根据偏导数的定义我们知道：
$\frac{\partial (\sigma^TE_i^Tu)}{\partial x}=\sigma ^T\frac{\partial E_i^Tu}{\partial x}+E_i^Tu\frac{\partial \sigma ^T}{\partial x}$

因此上式为

$(\frac{\partial \sigma^T}{\partial x}E_1^T+\frac{\partial \sigma^T}{\partial y}E_2^T+\frac{\partial \sigma^T}{\partial z}E_3^T)u\\=\frac{\partial \sigma^TE_1^Tu}{\partial x}+\frac{\partial \sigma^TE_2^Tu}{\partial y}+\frac{\partial \sigma^T E_3^Tu}{\partial z}-\sigma^T(\frac{\partial E_1^Tu}{\partial x}+\frac{\partial E_2^Tu}{\partial y}+\frac{\partial E_3^Tu}{\partial z})$

那么体力在整个体积上的积分则是：

$\iiint_\Omega [\frac{\partial \sigma^TE_1^Tu}{\partial x}+\frac{\partial \sigma^TE_2^Tu}{\partial y}+\frac{\partial \sigma^T E_3^Tu}{\partial z}-\sigma^T(\frac{\partial E_1^Tu}{\partial x}+\frac{\partial E_2^Tu}{\partial y}+\frac{\partial E_3^Tu}{\partial z}))]d\Omega$

显然我们可以看出来上式的第二项就是在整个体积上的应变能，因此有：
$\iiint _\Omega \sigma^T(\frac{\partial E_1^Tu}{\partial x}+\frac{\partial E_2^Tu}{\partial y}+\frac{\partial E_3^Tu}{\partial z})) d\Omega+(\frac{\partial \sigma^T}{\partial x}E_1^T+\frac{\partial \sigma^T}{\partial y}E_2^T+\frac{\partial \sigma^T}{\partial z}E_3^T)u \\=\iiint_\Omega (\frac{\partial (\sigma^TE_1^Tu)}{\partial x}+\frac{\partial (\sigma^TE_2^Tu)}{\partial y}+\frac{\partial (\sigma^TE_3^Tu)}{\partial z}) d\Omega$

也就是

$\iiint _{\Omega}\sigma^T E^T(\nabla)ud\Omega+\iiint_{\Omega}(E(\nabla)\sigma)^T ud\Omega\\=\iiint_\Omega (\frac{\partial (\sigma^TE_1^Tu)}{\partial x}+\frac{\partial (\sigma^TE_2^Tu)}{\partial y}+\frac{\partial (\sigma^TE_3^Tu)}{\partial z}) d\Omega$

最后根据高斯散度定理我们可以得到

$\iiint_\Omega \left\{\frac{\partial}{\partial x}(\sigma^T E_1^T u) + \frac{\partial}{\partial y}(\sigma^T E_2^T u) + \frac{\partial}{\partial z}(\sigma^T E_3^T u)\right\} d\Omega \\ = \iint_{\partial\Omega} \left\{n_x(\sigma^T E_1^T u) + n_y(\sigma^T E_2^T u) + n_z(\sigma^T E_3^T u)\right\} dS \\ = \iint_{\partial\Omega} \sigma^T(n_x E_1^T + n_y E_2^T + n_z E_3^T)u \, dS \\ = \iint_{\partial\Omega} [E(n)\sigma]^T u \, dS$

因此我们可以证明

$\iiint _{\Omega}\sigma^T E^T(\nabla)ud\Omega=\iint_{\partial \Omega}(E(n)\sigma)^T udS-\iiint_{\Omega}(E(\nabla)\sigma)^T ud\Omega$

这里可以看到为了推导出这个等式我们其实只用到了变形协调关系，平衡方程和边界条件，也就说只要满足了变形协调关系、平衡方程和边界条件，对一个连续的变形体这个等式就一定成立。注意，这里并没有用到本构关系，也就是说应力应变并不要求有对应关系，即力和位移也不需要有对应关系。我们知道在真实状态下，由于有外力所以产生了对应的位移和应力应变，但对于整个等式并不需要这个对应关系成立，只要变形协调条件、平衡方程以及边界条件成立即可。

到这里如果学过结构力学的单位荷载法我们就可以知道出单位荷载法为啥会成立了。这个方法主要是为了求变形体的位移。即如果要求一个变形体某点某方向上的位移，就在该点施加一个对应方向的单位荷载，称为虚拟状态，然后这个虚拟状态的单位荷载和实际位移相乘作为体力产生的外力势能就等于虚拟状态的应力和实际状态的应变相乘得到的应变能，因为是单位荷载，大小是1，可以不考虑，也就是说实际位移等于该虚拟状态下的应力和实际状态的应变相乘的应变能。

这个单位荷载法之所以成立就是对应了我上面说到的力和位移不存在关系，只需要满足分别满足变形协调、平衡条件和边界条件即可。因此我们考虑一个变形体受到的外力就是单位荷载，但有了一个实际位移。那么首先满足平衡条件和力的边界条件，对此我们只需要确定单位荷载下结构对应的应力即可，然后确定变形协调方程和位移边界条件，也就可以确定实际位移下对应的结构应变。因此单位荷载法就成立了。

那么根据这个等式我们可以推导出虚功原理

$\iiint _{\Omega}(\sigma^s)^T \delta \epsilon d\Omega=\iint_{\partial \Omega}\bar{p}^T \delta udS+\iiint_{\Omega}f^T \delta ud\Omega$

你可能已经发现这里就是把上面对应位置的位移和应变替换成了对应的变分形式，这里位移的变分形式也就称为虚位移，对应产生的应变也就是虚应变。

为啥要这么做呢？主要是为了推导后面的最小势能原理，因为最小势能原理的中的结论就是弹性体的位移函数使得总势能最小，那么就需要对这个位移函数进行变分。

3. 最小势能原理

下面我们就根据虚功原理来推导最小势能原理。首先定义总势能的表达式是：

$\Pi=\Pi_1+\Pi_2$

$\Pi_1 = \iiint\limits_{\Omega} Ud\Omega=\iiint\limits_{\Omega} \sigma^T E(\nabla)ud\Omega$

$\quad \Pi_2 = -\iiint\limits_{\Omega} f^T u d{\Omega} - \iint\limits_{S} \bar{p}^T u dS$

总势能是由弹性应变能和外力势能组成，这里可以看到外力势能是外力做的功的负值。那么有最小势能原理：在所有可能的位移中，弹性力学的精确解应该使得这个总势能达到最小。证明如下：

首先我们知道有实际位移状态 $u$ ，然后再定义一个可能的位移状态 $u^k$ ，正如我们上面说的这个可能的位移状态只需要满足变形协调和位移边界条件即可，但是这里进一步需要说明的是这个可能位移状态是一个任意接近于实际位移状态的函数，即：

$u^k=u+\delta u$

这就是变分的概念，关注的并不是自变量而是这里的位移函数本身。我们需要证明这个实际状态下的位移是可以让总势能函数达到最小的。对于线弹性体有应变能关系：

$U(u)=\frac{1}{2}\epsilon^T A \epsilon$

$U(u^k)=\frac{1}{2}(\epsilon^k)^T A \epsilon^k$

因此有：

$U(u+\delta u)=\frac{1}{2}(\epsilon+\delta \epsilon)^T A (\epsilon+\delta \epsilon)\\=\frac{1}{2}\epsilon^TA\epsilon+\epsilon^TA\delta \epsilon+\frac{1}{2}(\delta \epsilon)^TA\delta \epsilon$

$\Pi(u^k) = \Pi(u) + \iiint\limits_{\Omega} \varepsilon^T A\delta\varepsilon d\Omega + \iiint\limits_{\Omega} \frac{1}{2} (\delta\varepsilon)^T A\delta\varepsilon d\Omega - \iiint\limits_{\Omega} f^T \delta u d\Omega - \iint\limits_{\bar{\Omega}} \bar{p}^T \delta u dS \\ = \Pi(u) + \iiint\limits_{\Omega} \frac{1}{2} (\delta\varepsilon)^T A\delta\varepsilon d\Omega + \iiint\limits_{\Omega} (\sigma^T \delta\varepsilon - f^T \delta u) d\Omega - \iint\limits_{\bar{\Omega}} \bar{p}^T \delta u dS$