有限元 | 基于虚功原理推导梁单元刚度矩阵

fem178

发布于 2024-03-13 12:40:37

2260

发布于 2024-03-13 12:40:37

文章被收录于专栏：数值分析与有限元编程数值分析与有限元编程

梁的基本方程

▲图1

图1为受分布载荷作用的简支梁，该问题的三大类基本方程如下。

方向的平衡方程

-EI\frac {d^4 v}{d x^4} + p(x) = 0

方向的平衡方程

\begin{split} M(x) &= \int_A \sigma_x ydA \\ &= \int_A -y^2 E \frac {d^2 v}{d x^2} dA \\ & = -EI\frac {d^2 v}{d x^2} \end{split}

物理方程

\sigma_x = -Ey\frac {d^2 v}{d x^2}

几何方程

\epsilon_x = -y\frac {d^2 v}{d x^2}

单元的离散化描述

▲图2

图2所示为一局部坐标系中的2节点梁单元，其长度为

，弹性模量为

。单元节点位移列阵

\mathbf q^e = \begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \\ \end{Bmatrix}

其中

v_1,\theta_1,v_2,\theta_2

分别为各节点的挠度和转角。利用函数插值、几何方程、物理方程以及势能计算公式，可以将单元的所有力学参量用节点位移列阵

\mathbf q^e

及相关的插值函数来表示。例如，单元内部任意位置

处的位移(挠度)

可表达为:

v = \mathbf N \mathbf q^e =[N_1,N_2,N_3,N_4]\begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \\ \end{Bmatrix} \quad \cdots (1)

式中，

N_1,N_2,N_3,N_4

为Hermite插值基函数，其表达式为

\begin{split} N_1 & = 1-3(\frac {x}{l})^2 + 2(\frac {x}{l})^3 \\ N_2 & = x(1 - \frac {x}{l})^2 \\ N_3 & = 3(\frac {x}{l})^2 - 2(\frac {x}{l})^3 \\ N_4 & = x(\frac {x}{l}-1) \frac {x}{l} \\ \end{split}

▲图3

以单元节点位移作为独立自变函数，通过Hermite插值构造单元位移场。作为独立自变函数的位移首先要满足几何方程，位移边界条件以及单元间的连续性条件，故这种单元称为位移协调元。

单元应变场和应力场

由梁的几何方程，有梁的应变表达式

\begin{split} \epsilon &= -y \frac {d^2 v}{d x^2} \\ &= [\frac {d^2 N_1}{d x^2},\frac {d^2 N_2}{d x^2},\frac {d^2 N_3}{d x^2},\frac {d^2 N_4}{d x^2}] \mathbf q^e \\ &= -y[\frac {12x}{l^3}-\frac {6}{l^2},\frac {6x}{l^2}-\frac {4}{l}， \frac {6}{l^2}- \frac {12x}{l^3}，\frac {6x}{l^2}-\frac {2}{l} ] \mathbf q^e \\ &= \mathbf B \mathbf q^e \end{split}

其中，

\mathbf B

是应变矩阵。

由梁的物理方程

\sigma_x = E\epsilon = E \mathbf B \mathbf q^e

至此,单元上任意点的位移、应变和应力均已通过单元两端的结点位移表达。以下就利用虚功原理来导出单元的刚度矩阵。

虚功方程

设单元轴线处发生虚位移

\delta v

,则

\delta v = \mathbf N \delta \mathbf q^e

式中

\delta \mathbf q^e

为结点虚位移向量。则单元的虚应变

\delta \epsilon

可表示为

\delta \epsilon = \mathbf B \delta \mathbf q^e

存在于单元中的应力

\sigma

在上述虚应变中所作的虚功为

\begin{split} \delta W_i & = \int_V \delta \epsilon^T \sigma dV \\ & = \int_V \delta \mathbf q^{eT} \mathbf B^T E \mathbf B \mathbf q^e dV \\ & = \delta \mathbf q^{eT} \int_V \mathbf B^T E \mathbf B dV \mathbf q^e \end{split}

均布荷载

p_x

在虚位移上所作的虚功

\begin{split} \delta W_e & = \int_l p_x \delta v dx \\ & = \delta \mathbf q^{eT} \int_l p_x \mathbf N dx \\ \end{split}

由虚功原理

\delta W_e = \delta W_i

可得

\int_V \mathbf B^T E \mathbf B dV \mathbf q^e = \int_l p_x \mathbf N dx

记

\mathbf k^e = \int_V \mathbf B^T E \mathbf B dV

\mathbf F = \int_l p_x \mathbf N dx

则有

\mathbf k^e \mathbf q^e = \mathbf F

上式就是单元平衡方程。与最小势能原理相比，虚功原理不需做变分运算。

记

\mathbf B =[B_1,B_2,B_3,B_4]

，则

\begin{split} \mathbf k^e &= E\int_0^l \int_A y^2 \begin{bmatrix}B_1 \\ B_2 \\ B_3 \\ B_4 \end{bmatrix}\begin{bmatrix}B_1 & B_2 & B_3 & B_4\\ \end{bmatrix}dA dx\\ &= E\int_A y^2 dA\int_0^l \begin{bmatrix}B_1B_1 & B_1B_2 & B_1B_3 & B_1B_4\\ B_2B_1 & B_2B_2 & B_2B_3 & B_2B_4\\ B_3B_1 & B_3B_2 & B_3B_3 & B_3B_4\\ B_4B_1 & B_4B_2 & B_4B_3 & B_4B_4\\ \end{bmatrix}dx\\ & = \frac {EI}{l^3} \begin{bmatrix}12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\ -12 & 6l & 12 & -6l \\ 6l & 2l^2 & -6l & 4l^2 \\ \end{bmatrix} \end{split}

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2024-03-10，如有侵权请联系 cloudcommunity@tencent.com 删除

split