社区首页 >专栏 >有限元 | 弹性支座

有限元 | 弹性支座

fem178

发布于 2024-04-10 02:24:50

1500

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

考虑一段一端由线性弹簧支撑的细长梁，如图1所示，弹簧刚度为

，则弹簧的应变能为

▲图1

U_k = \frac 12 r(\omega_l)^2

其中，

\omega_l

是梁在

x=l

处的挠度。

梁的势能为

\Pi = \frac 12 \int_V \sigma \epsilon dV + \frac 12 r(\omega_l)^2 \quad \cdots (1)

现在来求

\omega_l

▲图2

如图2所示的两节点梁单元，

\omega_1,\theta_1,\omega_2,\theta_2

分别为四个节点自由度。物理坐标系

(xOy)

和自然坐标系

(\xi)

的线性映射关系

x = (1-\xi)x_1 +\xi x_2

在节点1的位置时

\xi=0

,在节点2的位置时

\xi=1

。挠度由三次多项式表示为

\omega = b_1 + b_2\xi + b_3\xi^2 + b_4\xi^3 \quad \cdots (2)

梁的转角

\begin{split} \theta &= -\frac {\partial \omega }{\partial x} \\ &= -\frac {\partial \omega }{\partial \xi}\frac {\partial \xi }{\partial x} \\ &= -\frac 1l(b_2+2b_3\xi + 3b_4\xi^2) \\ \end{split}\quad \cdots (3)

其中，

l= x_2-x_1

为单元长度。在节点1的位置时

\xi=0

,由(2)(3)可得

\omega_1 = b_1\quad \cdots (4.1)

\theta_1 = -\frac 1l b_2 \quad \cdots (4.2)

在节点2的位置时

\xi=1

,由(2)(3)可得

\omega_2 = b_1 + b_2 + b_3 + b_4 \quad \cdots (4.3)

\theta_2 = -\frac 1l (b_2 + 2b_3 + 3b_4 )\quad \cdots (4.4)

(4.1)~(4.4)联立解得

\begin{split} b_1 &= \omega_1 \\ b_2 &= -l \theta_1 \\ b_3 &= -3\omega_1 + 3\omega_2 + 2l\theta_1 + l\theta_2 \\ b_4 &= 2 \omega_1 -2\omega_2 -l\theta_1 - l\theta_2 \end{split} \quad \cdots (5)

(5)代入(2)得

\omega = N_1\omega_1 + N_2l \theta_1 + N_3\omega_2 + N_4l \theta_2 \quad \cdots (6)

其中

\begin{split} N_1 &= 1-3\xi + 2\xi^3 \\ N_2 &= -\xi + 2\xi^2 -\xi^3 \\ N_3 &= 3\xi^2 -2\xi^3 \\ N_4 &= \xi^2 -\xi^3 \end{split} \quad \cdots (7)

记

\mathbf N=[N_1, N_2l, N_3, N_4l]

,(6)写成矩阵形式

\omega = \mathbf N \mathbf q^e \quad \cdots (8)

由(3)可知，梁的广义应变(曲率)

\begin{split} \epsilon &= \frac {\partial^2 \omega }{\partial x^2} \\ &= \frac {\partial }{\partial \xi }(\frac {\partial \omega }{\partial\xi } \frac {\partial \xi }{\partial x }) \\ &= \frac {1}{l^2}\frac {\partial^2 \omega }{\partial \xi^2} \\ &= \frac {1}{l^2}(\frac {\partial^2 N_1 }{\partial \xi^2}\omega_1+ \frac {\partial^2 N_2 }{\partial \xi^2}l \theta_1 + \frac {\partial^2 N_3 }{\partial \xi^2}\omega_2 + \frac {\partial^2 N_4 }{\partial \xi^2}l \theta_2) \\ &= \frac {1}{l^2}[-6+12\xi, (4-6\xi)l, 6-12\xi, (2-6\xi)l]\begin{Bmatrix} v_1 \\ \theta_1 \\ v_2 \\ \theta_2 \\ \end{Bmatrix} \\ \end{split}

记

\mathbf B = \frac {1}{l^2}[-6+12\xi, (4-6\xi)l, 6-12\xi, (2-6\xi) ]

则

\epsilon = \mathbf B \mathbf q^e

其中

\mathbf B

是应变矩阵。应力

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

由(7)(8)可知，(1)中的

\omega_l

就是

\xi=1

时

\omega

的值。

\begin{split} N_1(1) &= 1-3 + 2 = 0 \\ N_2(2) &= -1 + 2 -1 =0 \\ N_3(3) &= 3-2=1 \\ N_4(4) &= 1-1=0 \end{split}

\omega_l = [0,0,1,0] \mathbf q^e \quad \cdots (10)

梁单元势能的表达式

\begin{split} \Pi &= U + U_r \\ &= \frac{1}{2} \mathbf q^{eT} \mathbf k_1^e \mathbf q^{e} + \frac{1}{2} \mathbf q^{eT} \mathbf k_2^e \mathbf q^{e} \end{split} \quad \cdots (11)

其中记

\mathbf k_1^e = EI \int_0^1\mathbf B^T \mathbf B d\xi

\mathbf k_2^e = r \mathbf N(\xi=1)^T \mathbf N(\xi=1)

\mathbf k^e = \mathbf k_1^e + \mathbf k_2^e

经计算得

\mathbf k^e = \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}+ \begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & r & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}

对于具有弹性支承的单元，只需将

\mathbf k_2^e

加入到常规单元刚度矩阵中。

▲图3

对于图3所示得梁，弹簧刚度为

r= \frac {EI}{l^3}

。只划分一个单元时，其有限元平衡方程为

\frac {EI}{l^3} \begin{bmatrix}12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\ -12 & 6l & 12+1 & -6l \\ 6l & 2l^2 & -6l & 4l^2 \\ \end{bmatrix} \begin{Bmatrix}\omega_1 \\ \theta_1 \\ \omega_2 \\ \theta_2 \\ \end{Bmatrix}= \begin{Bmatrix}R_1 \\ R_2 \\ -P \\ 0 \\ \end{Bmatrix}

考虑边界条件之后

\frac {EI}{l^3} \begin{bmatrix} 13 & -6l \\ -6l & 4l^2 \\ \end{bmatrix} \begin{Bmatrix} \omega_2 \\ \theta_2 \\ \end{Bmatrix}= \begin{Bmatrix} -P \\ 0 \\ \end{Bmatrix}

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

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

int